{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<h1><center>基于 Numpy 的面向数组有限元编程技术</center></h1>\n",
    "$$\\quad$$\n",
    "<h2><center>报告人: 魏华祎</center></h2>\n",
    "$$\\quad$$\n",
    "<h2><center>湘潭大学数学与计算科学学院</center></h2>\n",
    "$$\\quad$$\n",
    "<h2><center>weihuayi@xtu.edu.cn</center></h2>\n",
    "$$\\quad$$\n",
    "<h2><center>2017.04.08</center></h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 面向数组编程（Array-oriented programming)\n",
    "\n",
    "面向数组编程的核心思想是基于**多维数组及其运算**来实现相关算法。\n",
    "\n",
    "\n",
    "Users of computers and programming languages are often concerned primarily with the efficiency of execution of algorithms, and might, therefore, summarily dismiss many of\n",
    "the algorithms presented here. Such dismissal would be short-sighted, since a clear statement of an algorithm can usually be used as a basis from which one may easily derive more efficient algorithm. --- **Kenneth E. Iverson (Notation as a tool of\n",
    "    thought, 1979 ACM Turing Award Lecture)**\n",
    "\n",
    "**注:**\n",
    "\n",
    "肯尼斯·艾佛森（Kenneth E. Iverson，1920年12月17日-2004年10月19日）是一位计算机科学家，最重要的贡献是开发了APL。1979年他因对数学表达式和编程语言理论的贡献而得到图灵奖。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 面向数组编程的优缺点\n",
    "\n",
    "* **优点**:\n",
    "    - 避免写冗长的循环, 编写代码速度更高\n",
    "    - 代码简洁, 可读性更高\n",
    "    - 编程的抽象层次更接近高等数学的运算过程\n",
    "    - 可直接利用高效实现的数值代数软件包, 编写高效的代码\n",
    "* **缺点**: \n",
    "    - 需要改变已有的标量化编程思维方式\n",
    "    - 需要基于数组及其运算重构已有算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 支持面向数组编程的语言\n",
    "\n",
    "面向数组编程的语言或软件包:\n",
    "* APL\n",
    "* J\n",
    "* Fortran 90\n",
    "* Matlab \n",
    "* Octave\n",
    "* R\n",
    "* Wolfram Language\n",
    "* Python 中的 Numpy 扩展模块\n",
    "* C++ 中基于算子重载的很多数值代数扩展包, 如 Armadillo 和 Blitz++\n",
    "* ......\n",
    "\n",
    "\n",
    "下面分别是以 C 和 Python 语言实现两个矩阵的相加的代码:\n",
    "\n",
    "\n",
    "```C\n",
    "// C 语言实现 两个矩阵相加的代码\n",
    "for (i = 0; i < n; i++)\n",
    "    for (j = 0; j < n; j++)\n",
    "        a[i][j] += b[i][j];\n",
    "```\n",
    "\n",
    "```Python\n",
    "# Python 中两个矩阵相加的代码\n",
    "A += B\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  iFEM\n",
    "\n",
    "加州大学的陈龙教授开发了 Matlab 面向数组的有限元软件包 iFEM (https://bitbucket.org/ifem/ifem/), 利用很少的代码, 可以写出高效有限元代码. 但 iFEM 的主要缺点是**面向过程**编写, 导致代码的重用性和扩展性不高, 使用起来还是比较困难的.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FEALPy: Finite  Element Analysis Library in Python\n",
    "\n",
    "为了克服 iFEM 中存在的问题, 我们最近开发了 Python 语言的有限元软件包 FEALPy（https://github.com/weihuayi/fealpy）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "选择 Python 的原因如下:\n",
    "\n",
    "* 解释性语言，无需编译，易学易用。\n",
    "* 精心设计的语言, 允许编写结构良好且可读性强的代码(we code what we think)。\n",
    "* 支持面向对象软件开发模式，可用面向对象的方式组织有限元软件包，易重用、维护和扩展。\n",
    "* 有庞大的用户社群，包括 NASA, ANL, Google 等国际知名的科研机构和公司都把 Python 做为高性能计算的主要开发语言之一。Python 的初学者和开发者很容易从社群中获得帮助和开发文档。\n",
    "* 有丰富的科学计算基础软件包， 如：\n",
    "    + NumPy:  http://numpy.scipy.org - Numerical Python，主要提供多维数组及相关运算功能。\n",
    "    + SciPy:  http://www.scipy.org - Scientific Python，提供高效的优化、FFT,、稀疏矩阵等科学计算模块。\n",
    "    + Matplotlib: http://www.matplotlib.org - Graphics library，提供成熟 2D 和 3D 画图软件功能。\n",
    "* 对高性能计算有很好的支持，如MPI、OpenMP和CUDA等， 可直接在大型并行计算机上布署使用。\n",
    "* 很容易用其它编译语言扩展Python，如 C/C++和Fortran。\n",
    "* 开源免费。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "但 Python 也有一些问题, 如:\n",
    "\n",
    "* 开发环境不如 Matlab\n",
    "* 没有实现所有在专业软件包中找到的算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## FEALPy 的架构设计\n",
    "\n",
    "![FEALPy 的架构设计](./figures/sa-crop.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FEALPy 中已经实现的功能\n",
    "\n",
    "* 基础的网格数据结构、算法和可视化功能，包括二维的三角形、四边形和多边形，三维的四面体、六面体和多面体网格的数据结构；\n",
    "* 简单区域上的三角形网格生成算法，二维的界面拟合网格生成算法，曲面三角形网格生成；\n",
    "* 三角形网格的二分法加密算法。\n",
    "* 基于四叉树和八叉树的数据结构和加密粗化算法，在加密和粗化时，不再要求相邻叶子单元的层数之差最多为一层。\n",
    "* 实现了三角形和四面体网格上的任意次Lagrange有限元空间;\n",
    "* 多边形网格上的最低次虚单元空间;\n",
    "* Laplace 方程的离散、边界条件处理与求解等功能模块"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 面向数组的有限元编程技术\n",
    "\n",
    "有限元算法中的基本要素:\n",
    "\n",
    "* 网格\n",
    "* 有限元空间\n",
    "    + 基函数的计算\n",
    "    + 自由度的管理 \n",
    "* 离散系统\n",
    "    + 矩阵和向量的构造\n",
    "    + 求解\n",
    "* 可视化\n",
    "* ......"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 网格\n",
    "\n",
    "**目标:** 用 Numpy 中的数组来组织三角形网格中的数组结构.\n",
    "\n",
    "网格是有限元算法的基础, 最少需要两个二维数组来存储网格的信息:\n",
    "\n",
    "* 网格节点坐标数组 `point`:\n",
    "    + $N\\times 2$ 的二维数组\n",
    "    + `point[i, 0]` 和 `point[i, 1]` 分别存储第 $i$ 个网格节点的 $x$ 和 $y$ 坐标\n",
    "* 单元顶点编号数组 `cell`:\n",
    "    + $NC\\times 3$ 的二维数组\n",
    "    + `cell[i, 0]`, `cell[i, 1]` 和 `cell[i, 2]` 分别存储第 $i$ 个单元三个顶点的全局编号(即 `point` 中的行号)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "node = np.array(\n",
    "    [(0.0, 0.0),\n",
    "     (1.0, 0.0),\n",
    "     (1.0, 1.0),\n",
    "     (0.0, 1.0)], dtype=np.float)\n",
    "cell = np.array([\n",
    "        (1, 2, 0), \n",
    "        (3, 0, 2)], dtype=np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAEaCAYAAABARRODAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF5dJREFUeJzt3Xt01PWZx/HPk4QQ7hC5SqxY0kVAlLootm6tsFJaWtNVqWitwoK71bV11dV66lq8LF7as+5665bapRaxlfZQXdgKuIrH3lzwcBQvjbWJa5RLQARD5RKSTL77x0xwCIkJMDO/5zfzfp3DOblM8nsmTnxnfjN5YiEEAQDgVVHUAwAA8FEIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1QAQBcI1QAANcIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1QAQBcI1QAANcIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1QAQBcI1QAANcIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1QAQBcI1QAANcIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1QAQBcI1QAANcIFQDANUIFAHCNUAEAXCNUAADXCBUAwDVCBQBwjVABAFwjVAAA1wgVAMA1QgUAcI1Q5YCZzTGz30U9BwDEEaECALhGqAAArhEqSWZWZ2Y3mNkrZrbHzBaZ2TAzW2VmH5jZM2Y2KHXZM8zseTNrMLOXzezstM8zx8z+L/Uxb5nZJe2O869m9n7qfV/I8dUEgFgiVB+6QNI0SX8h6VxJqyTdJGmIkl+nq81spKQnJS2QVC7pekm/NLMhZtZH0v2SvhBC6Cfp05I2pH3+yZLekDRY0vckLTIzy8UVA4A4I1QfeiCEsC2EsFnSbyWtCyG8FEJolPSEpE9K+pqklSGElSGE1hDC05LWS5qR+hytkk4ys14hhPoQwh/SPv/bIYQfhRASkhZLGiFpWK6uHADEFaH60La0l/d18HpfScdL+krqtF+DmTVI+itJI0IIeyTNknSFpHoze9LMTkz7HFvbXggh7E292DcL1wMA8orbUJnZo2ZWb2Z/NrM/mdnlUc8kaaOkJSGEgWn/+oQQ7pakEMJTIYRpSt5b+qOkH0U57Ecxs56px+LeTj2mtoHHzQB45DZUku6SNCqE0F9SlaQFZvaXEc/0qKRzzWy6mRWbWZmZnW1mFaknX3w59VjVfkm7lTwV6FWJkuH9rKQBkm6W9AszGxXhTABwCLehCiH8IYSwv+3V1L/REY6kEMJGSV9W8kkW25X8H/0NSn4diyRdJ2mLpJ1KBuDKaCbtWghhTwjh1hBCXerxtl9JektS1D8MAMBBLIQQ9QydMrP/kDRHUi9JL0k6K4SwO9Kh8pSZDZP0tqSJIYQ/Rj0PALRxHSpJMrNiSZ+SdLak74YQmqOdKP+YWQ8ln47/Zgjh61HPAwDp3J76axNCSIQQfiepQo5PpcWVmRVJWiKpSdI3Ih4HAA5REvUAh6FEET9GlW9Sv3C8SMnf55rBvVUAHrm8R2VmQ83sIjPrm3p23XRJF0taE/VseeYHksZKOjeEsC/qYQCgIy4fozKzIZKWSTpFyZi+Len+EILb30uKGzM7XlKdkk+lb0l719dDCD+NZCgA6IC/UJkFKfkLSA9NGNr8ownDWrr4CHTD37+6reTyV9/tUdz2hhDYMwggFtw+RtXcu4f6X//pHlXTRveIepZ8sPvf/ld7X31X/aIeBAAOk7vHqBKS9vfuoc2njlDt1BOiHicv1D5Xp39+cJ3eHD1Ie4tN8r0xAwAO4u4e1e3Flhjzgy8W1049QaHYXUdjp/a5Ov183nJ97tYpWnHJBG2f+YvGB9ZtviPquQCgu9yV4I7iotaaaaOJVAYciNQtZ+u0y05RKC7ScyP6JRTCgqhnA4DuogZ5qn2kACCuCFUeIlIA8gmhyjNECkC+IVR5hEgByEeEKk8QKQD5ilDlASIFIJ8RqpgjUgDyHaGKMSIFoBAQqpgiUgAKBaGKISIFoJAQqpghUgAKDaGKESIFoBC5256OjhVypHr16rW1sbFxWHcuW1ZW1trY2BirH8DiNnNZWdm2ffv2DY96DhQOQhUDhRwpSWpsbBzW3b9EbWZF7v5qdRfiNrOZdeuHBiBTYvNTXKEq9Ei1t3r1ao0ZM0aVlZW6++673c6xf/9+zZo1S5WVlZo8ebLq6upyP2SauM0LpCNUjhGpgyUSCV111VVatWqVqqur9dhjj6m6utrlHIsWLdKgQYNUW1ura6+9VjfeeGPO52wTt3mB9ghVmnU/fkk//Pyj+pdR9+qJa1ZHOguROtQLL7ygyspKffzjH1dpaakuuugiLV++3OUcy5cv1+zZsyVJM2fO1Jo1axTV6b24zQu0R6jS9BveV2f942R98qKTIp2DSHVs8+bNOu644w68XlFRoc2bN2v+/PlasWJF5HN0dpmSkhINGDBAO3bsyNmMnc0i+Z8XaI8nU6QZN+MTkqQtr2xTc/3uSGYgUofv9ttvj3oEAFnEPSpHiNRHGzlypDZu3Hjg9U2bNmnkyJEu50i/TEtLi3bt2qVjjjkmp3N2NIvkf16gvYIOVVlTgyq3PqMzax7QWW/cozNrHlDl1mdUkmjM+SxEqmunnXaaampq9NZbb6mpqUlLly5VVVWVyzmqqqq0ePFiSdKyZcs0depUmVnOZ5XiNy/QXsGe+ivf/ZbGbVkhC60qUqskqaS1SSN2vaoRDUF7mo/r4jNkDpHqnpKSEj344IOaPn26EomE5s6dq/Hjx2v+/PmaNGlSzqLVnTnmzZunSy+9VJWVlSovL9fSpUtzMls+zAu0Z96e2VPSs6TpO3XX9MjmMcqaGjSpbrGKQ0uH71+0SHp3u2nqD+eqsXRgNkfJeaSW/cOTe1594vW+WT9QBplZOIxf+I3ds9XiNnNqXu5uIWcK8tRfxc71stB6yNsTCampSWptlVpbg4bWv6BEy6GXyxTuSQFA1wry1N+wD14/cLov3ZIlUuo0vSTpmadf1Wev66sp13864zMQKQDonoIMVXFrU4dvnzMn+a9NkPSbMUQqamVlZa1m1q17/2VlZbF7EkDcZi4rK8veaYYsiMNS4yiOG+F1PeylxgUZqkRRqUo6iVX7y2UakTp8jY2N3V7aGrfHe6T4zdzdHxq8iMNS4yiOG+F1PeylxrG6wWXKtn5j1drFVW9Vkbb1G5vR4xIpADh8BRmqTeWTFLr4oTBYkTaVT8rYMYnU0Zs7d66GDh2qk06KdsVVV3OEEHT11VersrJSJ598sl588cUcT3gotqdHt3k/qtvtxo0bNWXKFI0bN07jx4/Xfffdl5PjZuPrXJChaiwdqOpjq5SwkkPuWbWqSAkrUfWxVRl7ajqRyow5c+Zo9epolwV3Z45Vq1appqZGNTU1euihh3TllVfmcLpDsT092s37Ud1uS0pKdM8996i6ulpr167V97///axf52x9nQsyVJK0s+8JWj9qtuoHTFBLUamCpJaiUtUPmKD1o2ZrZ98TMnIcIpU5Z511lsrLy6Meo8s5li9frssuu0xmpjPOOEMNDQ2qr6/P4YQHY3t6tJv3o7rdjhgxQqeeeqokqV+/fho7duwhy4gzLVtf54J8MkWbxtKBqh1+jmqHn5OVz0+ksm/hwoWSpCuuuCLiST7U2bbyESNGuJln3bp1nV4mfXv64MGDczprtnT2Ncj1VpOo1NXV6aWXXtLkyZOzepzu3NaOREGHKpuIVG54ChTipxA27+/evVsXXHCB7r33XvXv3z/qcY5IwZ76yyYiVdi8bHk/nHnyfXu6t/8mudLc3KwLLrhAl1xyic4///ysHy9bX2dClWFEClVVVXrkkUcUQtDatWs1YMCAyE77SWxPl/xs3s+lEILmzZunsWPH6rrrrsvJMbP1debUXwYRqey6+OKL9dxzz+m9995TRUWFbrvtNjU3N0vK7SnAruaYMWOGVq5cqcrKSvXu3VsPP/xwzmbrCNvTo92839HtZd68eVk7Xpvf//73WrJkiSZMmKCJEydKku68807NmDEja8fs7Ot8tApye3o2xCVSbE/3J24zx217ehxuP1EcN+Lreli3H079ZUBcIgUAcUSojhKRAoDs4jGqo0CkcoPt6b7EbXt6HG4/URw3wut62LcfQnWEiFTusD3dlxhuT3d/+ymwx6gO+/YTqxucF0QKAHKHUB0mIhWdqLZBtxfH7elxnDnTotpi3tjYqNNPP12nnHKKxo8fr1tuuSVnx871xvhsfY8SqsNApKIVxTbojsRte7oUz5kzLaot5j179tSzzz6rl19+WRs2bNDq1au1du3arB83io3x2foeJVTdRKSiF8U26I7EbXu6FM+ZMy2qLeZmpr59k7+62NzcrObm5pw8iSGKjfHZ+h4lVN2Q7UjtfX+fls5drjtG36d/P+0hvfL46xk/Rr5J3wa9cOHCA1vUPehse7pncZw5E3J120kkEpo4caKGDh2qadOmZX2LuRT9f9NMbmznWX9dyMU9qZU3rVFxjyJd/8qV2vrau/rZZU9o+PghGjomP/7EQqa13wbNBnUcqVzddoqLi7VhwwY1NDTovPPO02uvvRb5X6rOpkxvbOce1UfIRaSa9jaremWNpnzrTPXsU6rjJ1dozOdG6+VluX/sJQ5yvQ36SMRxU3ccZ46jgQMHasqUKTl5rCyq/6bZ+B4lVJ3I1WNSO97cqaLiIg0e/eG582Hjhmj7Gzuydsy4imIb9JHwtj29O+I4c1xs375dDQ0NkqR9+/bp6aef1oknnpj140axMT5b36Oc+utANiNV1tSgip3rNeyD11Xc2qQNb5SoT5/k2xtLByYv07+n9u9pyuhx80Fn26DfeecdSbk7jRO37elSPGfOtKi279fX12v27NlKJBJqbW3VhRdeqC996UtZO16bbG0y/yjZ2tjO9vR2shmp8t1vadyWFbLQqiIlt4jU1Ejf/Kb05FMlqj62Sjv7nqDnF65X3fMb9dVHzsvo8SW2p3sUt5nZnp55BbaZgu3pRyPb96TGbVmh4tByIFKSVFEhJRJS/cYWjduyQmVNDdpavV1DxuTPX1cFgKNBqFKy/ZhUxc71snDoLsZevaTPfEZ6+GGpcW9Ce5/9td54qlanzByX8RkAII54jEq5eeLEsA9eP+ieVLprrpG+9z1p5vlB/frX6ot3zeCp6WnisP36aMRt5hhuT99mZsO6edlu39YyKYrjRnhdtx3uxxR8qHL17L7i1s6fHNG/v7RgQfLlIOk3Y8ZmbY44isP266MRt5njtj193759w6OeAUcnVje4TMvlWqREUWlGL1eIolzu2V5Xyz7379+vWbNmqbKyUpMnT1ZdXV3uh0zDUlrEWcGGKte7+7b1G6vWLr7crSrStn7cm+pMVMs92+vOss9FixZp0KBBqq2t1bXXXqsbb7wx53OmYykt4qwgQxXFgtlN5ZMUujhjEqxIm8on5WSeOIpquWd73Vn2uXz5cs2ePVuSNHPmTK1ZsybS03sspUWcFVyootqC3lg6UNXHVilhJYfcs2pVkRKW/D2qtl/6Rcc6Wu45f/58rVixImczdGfZZ/plSkpKNGDAAO3Y4XfbSNQLTIGPUlBPpoj6T3Xs7HuC1o+afdBmikRRqbb1G6tN5ZOIVDd0tNzz9ttvj3osAFlUMKGKOlJtGksHqnb4Oaodfk5kM+SD9OWeud5C3Z1ln22XqaioUEtLi3bt2qVjjvH7S9wspYVnBXHqz0ukcHSiWu7ZXneWfVZVVWnx4sWSpGXLlmnq1Kmuf1eKpbT5w8yCmVWmXv6JmS2Ieqajlff3qIhU/uhsuef8+fM1adKkrG+GbtPZss/0OebNm6dLL71UlZWVKi8v19KlS3MyW2dYSos4y+ultETqUCyl9SduM8dtKW2hMbMg6RMhhFoz+4mkTSGEmyMe66jk7ak/IgUg7szsODN73My2m9kOM3sw9fa5Zva6mb1vZk+Z2fFRz5pNeRkqIgUg7sysWNKvJL0taZSkkZKWmtmXJd0k6XxJQyT9VtJjEY2ZE3kXKiIFIE+cLulYSTeEEPaEEBpDCL+TdIWku0IIr4cQWiTdKWliPt+ryqsnUxCp/BSH7ddHI24zH8n2axyR4yS9nYpRuuMl3Wdm96S9zZS8x/V2robLpbwJFZHKX2y/RoHaKOljZlbSLlYbJd0RQvhpRHPlXGx+ivsoRApAHnpBUr2ku82sj5mVmdmZkhZK+raZjZckMxtgZl+JctBsi32oiBSAfBRCSEg6V1KlpHckbZI0K4TwhKTvKvnEij9Lek3SFyIbNAdifeqPSAHIZyGEdyT9TQdvXyJpSScfY2kvz8nacDkU23tURAoACkMsQ0WkAKBwxC5URAr5zsyeM7PLo54D8MLtY1R739+nFf/0P3rz13XqXd5Lf/3tz6h3eS8ilSFm9g1JcyRNkPRYvpzLBpB/3IZq5U1rVNyjSNe/cqW2vvauHv3qLxVag6bfNoVIZcYWSQskTZfUK+JZAKBTLk/9Ne1tVvXKGk351pnq2adUzfta1NKU0PFnVBCpDAkhPB5C+C9Jfv8+esyYWZ2ZXW9mr5jZLjP7uZmVpd73d2ZWa2Y7zWyFmR2b9nHTzOyPqY95UMktA+mft6AWkALtuQzVjjd3qqi4SINHlx94TGrM9EoV9yiOejSgKxdK+rykEySdLGmOmU2VdFfqfSOUXHOzVJLMbLCkxyXdLGmwpDclndn2yQpxASnQnstQNe1tVs9+pQc9ceITU0Zp/56mqEcDunJ/CGFLCGGnpP+WNFHSJZJ+HEJ4MYSwX9K3JX3KzEZJmiHpDyGEZSGEZkn3Stqa9vkKbgEp0J7LUJX27qHGXfsPeuLE/g+a1LNPadSjAV1Jj8xeSX2V3IB9YFloCGG3kqdcR6betzHtfSH9dX24gLTBzBok7dSHC0iBguAuVDcnWos+9mSNEk0JnXnV6Qcek9pavV1DxhwT8XTxZolWnV3/QbHMvhP1LAVmi5LBkSSZWR9Jx0jarOQut+PS3mfprysZra+HEAam/esVQng+N6MD0XMXqu8kQvE371+nz/Xtoff+tENNe5v1zgub9cZTtTpl5riox4stS7Tqaxf/Unev31Im6VYzK0k90F8sqTi18NLts0Bj7jFJf2tmE82sp5Kn79aFEOokPSlpvJmdn/r6Xy0pfVt8wS0gBdqz5JkGR8yClDxPMllS2h++CWYdfwi69sUg/UyyfqnXTbpN0i3tLnZbCOHW3E6WP8ysTtLlIYRnUq/fKqkyhPA1M7tC0g2SBkl6XtIVIYRNqct9XtL9koYpub9tgqQlIYT/TL3/UknfUvJe2S5JT4cQ5ubwqgGRchsqSa2SblEIC6IcJ28kT/fdqrZ70WmLKwHAM3+hAgAgjbvHqAAASEeoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBrhAoA4BqhAgC4RqgAAK4RKgCAa4QKAOAaoQIAuEaoAACuESoAgGuECgDgGqECALhGqAAArhEqAIBr/w9O0ryj76FbpQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from fealpy.mesh.TriangleMesh import TriangleMesh\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "tmesh = TriangleMesh(point, cell)\n",
    "fig, axes = plt.subplots(1, 3)\n",
    "tmesh.add_plot(axes[0])\n",
    "tmesh.find_node(axes[0], showindex=True, markersize=25, fontsize=12)\n",
    "tmesh.find_cell(axes[0], showindex=True, markersize=100, fontsize=12)\n",
    "axes[0].set_title('mesh')\n",
    "\n",
    "for ax in axes.reshape(-1)[1:]:\n",
    "    ax.axis('tight')\n",
    "    ax.axis('off')\n",
    "axes[1].table(cellText=node, rowLabels=['0:', '1:', '2:', '3:'], loc='center')\n",
    "axes[1].set_title('node', y=0.3)\n",
    "axes[2].table(cellText=cell, rowLabels=['0:', '1:'], loc='center')\n",
    "axes[2].set_title('cell', y=0.35)\n",
    "plt.tight_layout(pad=1, w_pad=1, h_pad=1.0)\n",
    "plt.savefig('/home/why/meshdata.pdf')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "除了上述两个基本数组之外, 有限元方法中还需要更多的网格数据, 如:\n",
    "\n",
    "* 边数组 `edge`\n",
    "    + 二维 $NE\\times 2$ 数组 \n",
    "    + `edge[i, 0]` 和 `edge[i, 1]` 分别存储第 $i$ 条边的起点和终点的全局编号(即对应 `point` 数组中的行号)\n",
    "    + 如果第 $i$ 条边是边界边, 则规定从 `edge[i, 0]` 看向 `edge[i, 1]`, 网格离散区域一定在左手边\n",
    "* 边与单元的相邻关系数组 `edge2cell`\n",
    "    + 二维 $NE \\times 4 $ 的数组\n",
    "    + `edge2cell[i, 0]` 和 `edge2cell[i, 1]` 分别存储第 $i$ 条边左右两个单元的全局编号(即对应 `cell` 数组中的行号)\n",
    "    + `edge2cell[i, 2]` 和 `edge2cell[i, 3]` 分别存储第 $i$ 条边在左右两个单元中的局部编号\n",
    "    + 如果是边界边, 则\n",
    "        - `edge2cell[i, 0] = edge2cell[i, 1]` \n",
    "        - `edge2cell[i, 2] = edge2cell[i, 3]`\n",
    "        \n",
    "**注: `edge` 和 `edge2cell` 可以从 `cell` 中构造出来.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  `edge` 和 `edge2cell` 和构造算法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbEAAAEXCAYAAAAjlXpCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuUVOWZ7/HvAy00CoaQyK0bQYRA00K3IBhnjJFExLuDZAh4JWgmZsWJnkwySeZEVIZR40rmwAwOTKKZJDiBZJIYUJBIDvFGNByNkgCJSxFINyBglAjaTTdVz/ljV0PZ6Utdd9Wu+n3W2mt1V+3a77vfZ3c/tW/PNndHREQkinoUugMiIiKZUhITEZHIUhITEZHIUhITEZHIUhITEZHIUhITEZHIUhITEZHIUhITkVCY2XfNbGGGn73TzB7KdZ+kc2Z2vpk1ZvjZEWbmZlaR6361pyQmIseY2U4zuyDX83bw2fPNLG5mh9tN52SyvHKUbmLP9otAIt5N7eK1JNPl5Ures6SISCf2uHt1oTshabnc3X9R6E4k056YiABgZsuBU4FHEt+y/9HMrjCzrWZ20MyeMLOazuZNvP4/Zva6mf3ZzJ4ys9oM+3KamT1pZofMbD3wwXbvX29mu8zsT2Z2e/JeoZn1MLOvmNn2xPs/MrMBWQxNwZnZl81sd2I8XjazS4F/Aj6ZGP/NifmGmtlqM3vTzF41s08nXr+ok/k/ZWa/Tyz3NTP7TIb962lm3zCzN8zsNeDSdu+fltgeDpnZL8zs/uS9QjP7sJn9KrGdbTaz81Nu3N01adKkCXcH2AlckPj5Q8A7wDTgBOAfgVeBXu3nTfr8PKAf0BtYBLyU9N53gYWJn88HGrvox7PAvyaWcx5wCHgo8d444DBwLtAL+AbQmtTvW4HngOrE5/8TWFHosc0iJmOABmBo4vcRwOnAnW1jkjTvU8B/AJVAPXAA+FjivY7mvzSxLAM+CrwLTOwoRh3FO+m9m4E/AMOAAcAvAQcqkuL5jUS8zgXeTopnFfAn4BKCHatpid9PSWV8tCcmIp35JLDG3de7eyvBP6E+wF919gF3/467H3L3IwT/NOvM7H2dzD408c07eTrJzE4FJgO3u/sRd38KeCTpc58AHnH3Z9y9BZhP8A+zzc3A/3b3xqR+fCKMiwzyJEaQjMeZ2QnuvtPdt7efycyGAX8NfNndm939JeAB4PrOFuzua9x9uweeBB4HPtJFX37WLl6fTrw+C1jk7g3u/iZwT1K/2uI5391b3P0ZYHXSMq8F1rr7WnePu/t64HmCpNYtJTER6cxQYFfbL+4eJ9gjqOpo5sQhpXsTh/HeJvjmDu0OBSbZ4+79203vJNp9K/Fzm11JPw9N9KOtX+8SfHNvMxx4uO0fLfB7gkQwqPtVLj7u/ipwG0Ey3m9mK81saAezDgXedPdDSa/topN4AZjZxWb2XOLw40GCxNFZvAD+pl28vp3UdkPSfO3j9WYiTm2S5x0O/G1yciTYWxvSRT+OURITkWTJezR7CP7BAGBmRnC4aHcH8wJcDVwJXAC8j+CwFwSHqtKxF3i/mZ2U9Nqp7d4/dkGImfUBPpD0fgNwcbt/tpXuvpuIcvcfuPu5BPFw4Ov85fjvAQaYWb+k106lk3iZWW/gJwR72IPcvT+wlvTjBUFMhrVrN/m9AWZ2YtJryfM2AMvbxeskd783lYaVxEQk2T5gZOLnHwGXmtnHzewE4B+AI8CvOpgXgnNhRwj2ik4E7s6kA+6+i+Bw0l1m1svMzgUuT5rlx8DlZvZXZtaLYA8l+R/vMuBfzGw4gJmdYmZXZtKXYmBmY8zsY4mk0ww0AXGC8R9hZj0A3L2BIDb3mFmlmU0AbgTaLqB4z/wE56d6E5w3O2pmFwMXZtjNHwGfN7NqM3s/8JW2N5LieWcinufw3ng+RBDP6Ym9+UoLbsFI6cpVJTERSXYP8LXEIZ3LCc5X/DvwRuL3yxPnod4zr5l9Efg+wWGk3cA2gosrujLU/vI+sZmJ964GzgbeBO5ILBsAd98K/D2wkuBb/mFgP0ECBVhMcM7lcTM7lOjH2RmNRnHoDdxLEIPXgYHAV4H/Sbz/JzP7TeLnOQR7wHuAh4E7/Pgl8e+ZP3HY8fMECegtgjFPPlfVkUfaxevhxOvfBn4ObAZ+A/y03eeuAc4h+IKzEPghiXglku+VBFdPHiDYM/sSKeYnS1wdIiISSWbWFzgIjHb3HYXuj3TPzH4I/MHd78h2WdoTE5HIMbPLzezExHmzbwC/4/iFJFJkzGyymZ1uwT18FxHsef0sF8tWEhORKLqS4JDZHmA0MNt1WKmYDQaeIDj0+2/AZ939xVwsWIcT88DM5gI3Ja4mEhGRPInqzX8ikqY+ffq83tzcnNK9UpWVlfHm5ua8HqkJqY19TU1Ng/PZRi6lEyNQnEBJTKRsNDc3D0r1yIuZ9cj3UZqQ2ojUDc7pxAgUJyjTc2IWFAv9kpn91szeMbMHzWyQmT2WVKDy/Yl5Oy1MaWZzLSiaecjMdpjZNe3a+YaZvZV47+KQV1OkQ+vWrWPMmDGMGjWKe+9N6X7SjDQ0NDB16lTGjRtHbW0tixcvzks7Ya1P2MJcr3y3lddtIR8FK4t9IriK6TmCMjRVBPeY/AY4k6Bw5gaCe1M6LUwJnERQxHJMYplDgNrEz3MJCpJ+GugJfJbgBLQVet01le8E+NGjR33kyJG+fft2P3LkiE+YMMG3bt3q7QX/GrKzZ88ef+GFF9zd/e233/bRo0e/p61ctNHd+iTaKPjYpzq1jUmYcUpxDLOS4raQ0ZiV5Z5Ywr+7+z4PStE8Dfza3V9092aCmwTPpPvClHHgDDPr4+57PbgJs80ud/+2u8eA7xEkuUgd2pDSs2nTJkaNGsXIkSPp1asXs2fPZtWqVXlpa8iQIUycOBGAfv36UVNTw+7dua38FOb6hCnM9QqjrXxuC+WcxPYl/dzUwe996aIwpQfFST9JUDF7r5mtMbOxSct4ve0HP174sm8e1kMkZbt372bYsONl66qrq9m9ezfz589n9eruijVkbufOnbz44oucfXZuC2d0tj5RF2acwh7DXG8LJXlhhwUPW/s4wSG/14H73P2BDBbVVpjy0x296e4/B36eKEC6kKD0SlePMUhZok7afxAUUx0AbAe+6u6P5WL5IskWLFiQt2UfPnyYmTNnsmjRIk4++eS8tVMO8hmnMORjWyjVPbF7gBHufjJwBbDQzCZlsJxOC1MmLgS5MlEx4AjBTXzx3K0CFQRJ9KMEFcG/BvzIzEbksA0pM1VVVTQ0HH8KRmNjI1VVnT6pI2utra3MnDmTa665hquuuirnyw97fcIS5nqF1VbetoVMT6ZFZSJ4KupeYFbSaztJekIpQbK6M+n3m4BfJH4+G3iSoBDpAWANwWMGhiRe/zNB3bYngHGJz8wFnmnXDwdGZbkuvwVmFnpMNUVzAry1tdVPO+00f+21146dxN+yZYu3Rw5O5sfjcb/uuuv81ltv7fD9XLTR3foQ0Qs7woxTimOYlRS3hczGLNMPFvtEcCju3UTy+A3Qt9B9ynJ9BhE8hmFsofuiKZpT2z+jNWvW+OjRo33kyJG+cOFCd3e//fbbfdWqVd4mF/+4nn76aQd8/PjxXldX53V1db5mzZqctuHe8fq0a6PgY5/qlDwmYcWps7Zy2UaK20JGY1bSZafMrCdB+f/zga978Ij1yEk8y+kxYLu7f6bQ/ZFoMjNP9e/dzMj3/4YQ28jkIY8FkU6MEvOXfZxK9ZwYAO4ec/dnCJ4C+9lC9ycTiQfYLQdagFsK3B0RkaJSklcndqACOL3QnUhX4nHwDxIcSrwkqnuSIiL5UnJJzMwGAh8DHiW43+sCgqedzilkvzK0FKghuAilqdCdiTIzc4KHJr5qZt8FGt39awXuVqgqKyvjSY+m725egu9Qee1PGG3k8orhvEsnRon5yz5OJZfECC7k+CywjOBw6S7gNnfP352ceWBmw4HPEFy+/3rSRvQZd//vgnVMIqu5uTnlQq4ldK4lUqdM0okRKE5Qas8TC75tEzd49IY6HrtmQqF7lLKLH/otl31v8/GTlBE6GR0F2hM7ftHAvHnzePTRRxk4cCBbtmzpbN6c/ONat24dt956K7FYjJtuuomvfOUrOW+jq/WJ6oUdqcQoMX8k4tTd+mQTp5JMYkdOPIGfLL2UV6ZF5zRYxZ2/9M9/6zfWr+2FCP3hhc3MhgGLCaqj9ABWuPstZjYP+BLBU2Q3AX/n7rsSn1ESS/yDfOqpp+jbty/XX399XpNYLBbjQx/6EOvXr6e6uprJkyezYsUKxo0bl7M2gC7XJ6pJLJUYJeaPRJy6Wx9dnZgQJ0hguycO4dWPnVbo7qRs84+3+j9/f7PtrPkgTb17Qm4rf5SUxG0TjxIcJh5B8KSBlWZ2JfBPwFUETxl4GlhRoG4WtfPOO48BAwbkvZ2witiGtT5hCnOdwohTPtenpJLYt+sGHf3J0kt5aMVMvGc0Vm3zj7f6o1/+hc1YehkPP34d//6ZSRA8BkY6NgUYCnzJ3d9x9+bEbRQ3A/e4++/d/ShwN1CfOLco3Vi2bBnLli3L6TJLtThvoeQjRhD9OJXUhR0PTBh09LJpp0dmndoS2Mz/uJSx00cFpUXqBoP7wkL3rYgNI3jMzdF2rw8HFpvZN5NeM4I9tV1hdS6qbr755kJ3QbqhGHUsMv/wS037BCYpawBONbOKdomsAfgXXblZPEq1OG+piXqconHMrcQogWVlE0FB53vN7KTEkwX+muCWiq+aWS2Amb3PzP62kB0td5MnT+aVV15hx44dtLS0sHLlSq644opCd0vaiXqclMRCpgSWHQ+elH05MAr4I9AIfNLdHwa+TnCRx9vAFuDignW0iM2ZM4dzzjmHl19+merqah588MG8nG+pqKhgyZIlTJ8+nZqaGmbNmkVtbW1O24CO1yfqwooRhBOnfMaopC6xn3x9XdNlX59Wme1yfv2dF3npR1vZ/4c3OONvxjJj0UW56F5KCez3j73CynmrInNJsESHCgAXPxUATp/OiXWg3+C+nHfr2Wx/chetze2vH8iM9sBERHJPSawD4y4ZDcCe3+6jde/hrJenBCYikh9KYkBly0Gq33yeQYd+T894C7EevdjXr4YnY9kvuxwSWJ8+fV5vbm4elOr8lZWV8ebm5ozPx2b7+RwuY19TU9PgbJYRJhUALn4qAJy+sk9iAw7vYNye1ZjH6ZEolFERb2HIn3/HkIPOO63DullC58ohgQE0NzcPSvM4flpFTnP9+RwuI+XEXQxUALj4qQBw+iIV4FyrbDnIuD2r6elHjyWwNj2I0wOn/7sNVLYcTHvZ5ZLA2lu3bh1jxoxh1KhR3HvvvRkvZ968eQwcOJAzzjgj42U0NDQwdepUxo0bR21tLYsXL057Gblan2KSi3FJRS5i2J2w1iVsYa5X1ONU1kms+s3nMe9uL9apfvP5tJZbrgksFovxuc99jscee4xt27axYsUKtm3bltGy5s6dy7p167LqT0VFBd/85jfZtm0bzz33HPfff39a/cnl+hSTbMclVbmIYXfCWpewhbleUY9TWSexQYd+/xd7YACxGLS0QDwOHof3v7GN2NHUDtmWawKD3BYSzUXB0CFDhjBx4kQA+vXrR01NTVo14cIqYBu2bMclVWEUsQ1rXcIW5npFPU5lncR6xls6fH35cpg+HX7wA1i/Hi69sJWnFj3X7fLKOYFB54VE58+fz+rVhX0m6c6dO3nxxRc5++yzU/5M1AujpiJ5XPJ1M21YMolxFJRSjCD3cSrrCztiPXpR0UEimzs3mNoc7dGLjaP/qstllXsC68qCBQsK2v7hw4eZOXMmixYt4uSTTy5oX4pJ+3GJcoHZUo1xKcUI8hOnst4T29evhng3QxCnB/v61XQ5jxJYoBgLiba2tjJz5kyuueYarrrqqrQ+W4zrkyvZjEuxKaV1SVZq65Wv9SnrJNY44Cy8mys73XrQOOCsTt9XAjuu2AqJujs33ngjNTU1fOELX0j788W2PrmS7bgUk1Jal2Sltl75XJ+yTmLNvfqzbegVxKziL/bI4vQgZhVsG3oFzb36d/h5JbD36qyQaCbnxHJRMHTjxo0sX76cDRs2UF9fT319PWvXrk3582EVsA1bZ+OS6/MtYRTmzTbGxSqsGEH046QCwHResaNxwFmhJ7AoFgAOu2hpLm6+zOEyIhMrFQAufioAnL6yvrCjTXOv/rw6+AJeHXxBSvNrD0xEpDiU9eHETCiBiYgUDyWxNCiBiYgUFx1OTFG6Cezdt5pY/Q+Ps/3JnZw4oA8f/+pHmHBV15fqR1XYlbdzUVU7R8so2Qrpqo5eGKpinz4lsRRksge29p/+Lz1P6MEXf/tZXt+ynx9c/zCDa09h4JgP5rm34Qu78nYRXdgRqSMZqmJf/FTFPn2RCnAhZJLAWt5tZdvaV5j6j39N75N6MfzsasZceDqbfxz9wqRdyWU17ObmZqZMmUJdXR21tbXccccdGS0nmyr0pVohPVdjm4p8PwVAMcqNKMdJe2JdSDWBtb9E/+VXK6jo6VQP60FzYp5B405h17ON4XS8QObOncstt9zC9ddfn/WyevfuzYYNG+jbty+tra2ce+65XHzxxXz4wx9OeRltVejXr19PdXU1kydP5oorrmDcuHEpfb6t8vbEiRM5dOgQkyZNYtq0aSl/vljlYmxTke34p0Ixyl7U46Q9sU6kmsAGHN7BWTu/x5A//46KeAsGtLx7lJP6OGft/B4DDu8AoPLk3hx5p+OCw6Uil9WwzYy+ffsCQbma1tbWtI/LZ1uFvlQrpOdibFMRxlMAFKPsRT1OSmIdSGcPrKOHavbpA+++Cz39KOP2rKay5SBHDrXQ+6ReYXS/qGRTYSAWi1FfX8/AgQOZNm1a2lWvc1mFvtQqpHc0trl+2kDYTwFQjDIT9TjpcGI76ZwD6+yhmtXVwTPJGhthaHWc6jef5/VtrZwy5gP56nbRyqbqds+ePXnppZc4ePAgM2bMYMuWLXl9+mxnSrFCekdjW+inDWRDMYoGVbHPs3Qv4ujsoZp9+sBHPgL/9V9wpCnOgV9t5eWfv0rdJ6J9nL5Q+vfvz9SpU9N++mwuqtCXWiXx9jId21SE9RQAxSg7UY+TklhCJlchdvZQTYDbboMjR+Cqq+CeBUe59J4LSvLy+nw5cOAABw8eBKCpqYn169czduzYtJaRbRX6Uqsk3iYXY5uKMJ4CoBhlL+px0uFEMq/E0dlDNQFOPhkWLgx+Dh6qWZo3OiebM2cOTzzxBG+88QbV1dXcddddtLa2AukfVty7dy833HADsViMeDzOrFmzuOyyy9JaRnIV+lgsxrx589KqQt9WeXv8+PHU19cDcPfdd3PJJZek1Y9i09nYzp8/n7POOitn/8CyHf9UKEbZi3qcyr6KfTalpEa9/guG/Pl3HR5SbBOnB3vfNz7l4sKqYp//z+d4GZGJlarYFz9VsU9fWR9OzLYWYi4eqikiIpkr2ySWi2K+2T5UU0REslOW58RyWY3+zb6n8fyIG9J+qGYpqays3Gdmg9KYP60ip7n+fA6XsS+bz4dNBYCLnwoAp6/sklg+HqeS7kM1S01TU9PgQvdBuqcCwMVPBYDTF6kAZ0vPAxMJxGIxzjzzzLSv+ExHGMV5wy6UGzbFqXsluSfW0bO8PB6PXAIzs1uAucB4YIW7zy1oh6RkLF68mJqaGt5+++28tRFGcd4wC+UWguLUvZLcE0t+ltdVSy5h9Rcf55EvrY9UAkvYAywEvlPojkjpaGxsZM2aNdx00015bSeM4rxhFsoNm+KUmpJLYu2f5XWw4c8ea40x+uOnRS2B4e4/dfefAX8qdF+kdNx2223cd9999Ohx/M8/H4Vlk+WzOG+2haKLleKUmpJLYn/a/iY9evbgg6cPOHYOrG7mOGItkbpISSQvHn30UQYOHMikSZPe8/qCBQtyXmqoTb6L87YVym1sbGTTpk1s2bIl522ETXFKXcklsZZ3W+ndr9d7LuIYNnloyT/LSyQVGzduZPXq1YwYMYLZs2ezYcMGrr322ry1F2Zx3nwXyg2T4pS6kktivU48geaDzSRfxFGuz/ISae+ee+6hsbGRnTt3snLlSj72sY/x0EMP5aWtMIrzhlkoN0yKU+pKKondtHlfxaCVv/NYa5wL53/02Dmw17cdiMSzvCwWZ+Lm18Hs9kL3RcpLPs61tBV93bBhA/X19dTX17N27dqctrF3716mTp3KhAkTmDx5MtOmTcvr5eiFpjj9pZIqABwz83eBWSf3ovH807jiX6fz+pb9/Pd1P+XG1XOK+lEoFotz7ZyfMOT/7abPkVgc955mVkFwG8QdQDXwaeCoux8taGclklQAuPipAHAGny2lJIaZAzQYnHviCbGGpqM9e1SY969+X1PfU05sLXT3ujLtzaaKRdvfOumkeCIe7mZmdxIksGR3ufud4fZOSoGSWPFTEsvgs6WYxIA4cAfuCwvZnbQEhxDvpO0Qb4T+8CQalMSKn5JYBp8tqSQmIp3q06dPrLm5OeUCwM3NzXntT0htxJuamnrmtZEcSidGoDiBkphI2dCeWPHTnlj6SurqRBERKS9KYiJlaMSIEYwfP576+nrOOis/Tx4Pq8L8unXrGDNmDKNGjeLee+/NSxuFEEaM2uR7DPO6Lbi7Jk2aymAK/twDw4cP9wMHDnhnkufNVDwe90OHDrm7e0tLi0+ZMsWfffbZnLZx9OhRHzlypG/fvt2PHDniEyZM8K1bt7Zvo+Bjn+qUToyS1i8rKY5hVlLcFjIaM+2JiUhehFFhftOmTYwaNYqRI0fSq1cvZs+ezapVq3LaRqkLYwxVxV5EcsrMuPDCC5k0aRLf+ta3AFi2bBnLli3LaTv5rjC/e/duhg0bduz36urqnD9GpFDCilFYY5ivbaEkH4opIl175plnqKqqYv/+/UybNo2xY8dy880357ydtsrlBw8eZMaMGWzZsoUzzjgj5+2UorBiFJZ8bQvaExMpQ1VVVQAMHDiQGTNmsGnTpry2l68K81VVVTQ0NBz7vbGx8di6RV1YMQp7DFXFXkSy8s4773Do0KFjPz/++ON52TsKo8L85MmTeeWVV9ixYwctLS2sXLkyb8/bClNYMYJwxjCf24IOJ4qUmX379jFjxgwAjh49ytVXX81FF1107FxLrg5Z7d27lxtuuIFYLEY8HmfWrFk5rzBfUVHBkiVLmD59OrFYjHnz5lFbW5vTNgohrBhBOGOYz21BFTtEyoQqdhQ/VexInw4niohIZCmJiYhIZOmcmEiZqKys3Gdmg1KcN25mef2SG1Ib+/K5/FxLJ0aJ+cs+TjonJiIikaXDiSIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYiIiEllKYt0ws61mdn6h+yGdU4yiw8zON7PGQvcjysxsp5ldkKdlZxwfMxthZm5mFbnuV1eUxLrh7rXu/gSAmf0vM3vdzN42s++YWe8Cd084HiMzO8PMfm5mb5iZF7pfIlGXSJhNZnY4aVpS6H4lUxJLkZlNB74CfBwYDowE7ipop6S9VuBHwI2F7ohICbnc3fsmTbcUukPJlMS6kbTrfgPwoLtvdfe3gH8G5ha0cwIcj5G7v+zuDwJbC92ncmNmQ83sJ2Z2wMx2mNnnE6/3MbPvmtlbZrYNmNzucxPN7EUzO2Rm/2NmPzSzhUnvX2ZmL5nZQTP7lZlNCHnVipaZ9TazRWa2JzEtSj46ZGZXJsbubTPbbmYXJV7/lJn9PjHmr5nZZzJsv6eZfSNx5OM14NJ2759mZk8l2vmFmd1vZg8lvf/hREwPmtnmTE8JKImlrhbYnPT7ZmCQmX2gQP0RKQpm1gN4hOBvoorgaMVtiaMXdwCnJ6bpBF8G2z7XC3gY+C4wAFgBzEh6/0zgO8BngA8A/wms1mH8Y/438GGgHqgDpgBfAzCzKcD3gS8B/YHzgJ2Jz+0HLgNOBj4F/B8zm5hB+59OLOdM4CzgE+3e/wGwiSB2dwLXtb1hZlXAGmAhQey/CPzEzE5JtxNKYqnrC/w56fe2n/sVoC8ixWQycIq7L3D3Fnd/Dfg2MBuYBfyLu7/p7g3AvyV97sNABfBv7t7q7j8l+KfX5u+A/3T3X7t7zN2/BxxJfE7gGmCBu+939wMEpzfaEsWNwHfcfb27x919t7v/AcDd17j7dg88CTwOfKSLdn6W2Ftqmz6deH0WsMjdG9z9TeCetg+Y2akE28X8xDbxDLA6aZnXAmvdfW2if+uB54FL0h2EUK8iibjDBN9c2rT9fKgAfREpJsOBoWZ2MOm1nsDTwFCgIen1XUk/DwV2u3vyRTjJ8w4HbjCzv096rVficxKMQ/J47uL42AwD1nb0ITO7mGAP+UMEOzInAr/rop2/cfdfdNJ+V7F9093fTXqtIdEvCGL7t2Z2edL7JwC/7KIfHdKeWOq2Euyyt6kD9rn7nwrUH5Fi0QDscPf+SVM/d78E2Mvxf1wApyb9vBeoMjNLei153gaCvbjk5Z7o7ivytibRsocgGbQ5NfEaBGN3evsPJA7F/gT4BjDI3fsTJDtrP28KuovtADM7Mem19rFd3i62J7n7vel2Qkksdd8HbjSzcWbWn+DY83cL2yVJZoFKgm/rmFmlzp+EYhNwyMy+nLiQo2fidofJBFeLftXM3m9m1UDyXtWzQAy4xcwqzOxKgvM6bb4N3GxmZydie5KZXWpmOoQfWAF8zcxOMbMPAvOBtgsnHgQ+ZWYfN7MeZlbzdRzKAAANIklEQVRlZmMJ/jZ6AweAo4m9sgszbP9HwOfNrNrM3k9w9TYA7r6L4PDgnWbWy8zOAZL3uh4CLjez6YntpdKCe9Sq0+2EkliK3H0dcB/B7u4fCXad7yhop6S94UATx69ObAJeLlx3yoO7xwhO8NcDO4A3gAeA9xGcp9mVeP1xYHnS51qAqwjO3xwkOE/yKMF5L9z9eYKLB5YAbwGvoiuCky0kSBS/JTgc+JvEa7j7JhIXbRCcv38SGO7uh4DPEySgt4Cree+5qo480u4+sYcTr38b+DnBBT2/AX7a7nPXAOcAf0r064ccj20DcCXwTwQJtYHgIpS0c5K993C0iEjhmNmvgWXu/l+F7ovklpn9EPiDu+f0y7/2xESkYMzso2Y2OHE48QZgArCu0P2S7JnZZDM7PXE48yKCPa+f5bodXZ0oIoU0huDQ1knAa8An3H1vYbskOTKY4BDjB4BG4LPu/mKuG9HhRBERiSwdThQRkchSEhMRkchSEhMRkchSEhMRkchSEhMRkchSEhMRkcgqqfvE+vTp83pzc/OgQvcjW5WVlfuampoGF7of+aI4FUY6415ZWRlvbm7O65fckNqIVIwkfSV1n5iZeSmsj5nh7plUlY4Exakw0hn3xLrluz9htRGZGEn6dDhRREQiS0ksRPPmzWPgwIGcccYZhe6KdKIcYrRu3TrGjBnDqFGjuPfetB/fVHRtlUPMpHNKYiGaO3cu69aptmkxK/UYxWIxPve5z/HYY4+xbds2VqxYwbZt2yLdVqnHTLqmJBai8847jwEDBhS6G9KFUo/Rpk2bGDVqFCNHjqRXr17Mnj2bVatWRbqtUo+ZdE1JTKSM7N69m2HDjj8lvrq6mt27dzN//nxWr+7u2Yi5aUskl0rqEnsRycyCBQsK3QWRjGhPTKSMVFVV0dDQcOz3xsZGqqqqIt+WlC8lMZEyMnnyZF555RV27NhBS0sLK1eu5Iorroh8W1K+lMRCNGfOHM455xxefvllqqurefDBBwvdJWmn1GNUUVHBkiVLmD59OjU1NcyaNYva2tq8nBPrrK1cK/WYSddUsaMIlXqVAcWpMFSxQ0qR9sRERCSylMRERCSySuoS+8rKyriZRT4xV1ZWxgvdh3xSnAojnXGvrKzELL9H4UJqI1IxSvcJD3ragM6JFaVSP46vOBWGzokVv3T/NhQnHU4MXZjFVyV9zc3NTJkyhbq6Ompra7njjjsK3aWcC6tgbkNDA1OnTmXcuHHU1tayePHivLRTqn9TYa1XGNtDXrcFdy+ZKVid4nX06FEfOXKkb9++3Y8cOeITJkzwrVu3/sV8ifUo+HjmayrmOMXjcT906JC7u7e0tPiUKVP82Wef7XDeqMWpbdyffPJJf+GFF7y2trbTcchFjPbs2eMvvPCCu7u//fbbPnr06Pds77loo7u/qajGKM3/FVnpbnsIcVvIaMy0JxaiMIuvSmbMjL59+wLQ2tpKa2tr3s/bhC2sgrlDhgxh4sSJAPTr14+ampqc104s1b+pMNcrjO0hn9uCkliIVBA1GmKxGPX19QwcOJBp06Zx9tlnF7pLebds2TKWLVuWt+Xv3LmTF198MedjWap/U2EWag5brreFkro6USQXevbsyUsvvcTBgweZMWMGW7ZsKfkHLt588815W/bhw4eZOXMmixYt4uSTT85bO+Ug6oWa87EtaE8sRCqIGi39+/dn6tSpeuBiFlpbW5k5cybXXHMNV111Vc6XX6p/U6W4XvnaFpTEQqSCqMXvwIEDHDx4EICmpibWr1/P2LFjC9yraHJ3brzxRmpqavjCF76QlzZK9W+q1NYrr9tCpleEFONEEV/11mbNmjU+evRoHzlypC9cuLDDeYjYFVXpTsUcp82bN3t9fb2PHz/ea2tr/a677up03qjFqW3cZ8+e7YMHD/aKigqvqqryBx54wJcuXepLly5tv25Zefrppx3w8ePHe11dndfV1fmaNWty2oZ7139TUY1RZ+t1++23+6pVq9qvX1Y62h6ShbgtZDRmutm5CEXtBs10KU6FoZudi59udk6fDieKiEhkKYmJiEhkldQl9iosGw2KU2GoAHDxS/dvQ3FSAeCiFLXj+OlSnApD58SKn86JpS/y34ajZsSIEYwfP576+nrOOuusQndHOlDqMQqrMG9YxZRLsQBwWEWaoQQKNWd6WWMxThTxpdtthg8f7gcOHOhyHiJ2WXC6U7HHKZUYuUcvTm3j3l0x1qR1y0p3xZRz0UapFgBOpUhz0vplJeqFmrUnJlJmwijMC+EUUy7VAsBhFWmG6BdqVhILmZlx4YUXMmnSJL71rW8VujvSgXKKUXIx1nwUAc53MeVSLQDckXwXaYZoFmouqasTo+CZZ56hqqqK/fv3M23aNMaOHct5551X6G5JknKJUftirPkoAlyOxZTzJZ9FmiG6hZq1JxaytiKeAwcOZMaMGWzatKnAPZL2yiFG+S7M216+iimXYqHcQohyoWYlsRC98847HDp06NjPjz/+uL6VFplyiJF7/gvzQjjFlEutUG4hhLE95DNOOpwYon379jFjxgwAjh49ytVXX81FF11U4F5JsnKI0caNG1m+fPmx2wgA7r77bv74xz8CuTtstXfvXm644QZisRjxeJxZs2Zx2WWX5WTZbSoqKliyZAnTp08nFosxb948amtrc9pGIcyZM4cnnniCN954g+rqau666y5aW1uB3B9W7Gx7uOSSS3LWRj7jpJudi1DUbtBMl+JUGLrZufjpZuf06XCiiIhElpKYiIhElpKYiIhEVkld2KHq6NGgOBWGqtgXP1WxT58u7ChCUTsZnS7FqTB0YUfx04Ud6Yv8t+EoCatatGSuHGIUVnX5MCqxl2q8wl6vfD8JIK/bQqaVg4txosiro6dSPdw9u4rOUZiKOU6pxsg9enFqG/fuqssnrVtWuqvEnos2UqzAXvCxT3UijScNJK1fVlKsMJ+VFLeFjMZMe2IhCqt6uGSuHGIURnV5CKcSe6nGK8z1CuNJAPncFpTECiRf1aIld0o5Rh1Vl58/fz6rV68udNcyVqrxyveTBqL+JICSujoxKqJaLbqclHqMOqouv2DBgkJ3K2OlGq8wnjQQddoTC1nY1cMlfeUUo3xVlw9TqcYrrPWK+pMAlMRC5B5O9XDJXDnEKIzq8mEp1XiFuV6RfxJApleEFONEEV/15u7+9NNPO+Djx4/3uro6r6ur8zVr1vzFfETsiqp0p2KOU6oxco9enNrGffPmzV5fX+/jx4/32tpav+uuu9zd/fbbb/dVq1Ylr1tWZs+e7YMHD/aKigqvqqryBx54oKPxy0p38YpqjDpbr6VLl/rSpUtzOobu7mvWrPHRo0f7yJEjfeHChe95L8RtIaMx083ORShqN2imS3EqDN3sXPx0s3P6dDhRREQiS0lMREQiq6QusVdh2WhQnApDBYCLX2Vl5T4zG5TG/Hn/WwqpjX2ZflbnxIpQ1I7jp0txKgydE5NSFPlvw1ETi8U488wzueyyywrdFelCqccpjPULo4htWMWMpXgpiYVs8eLF1NTUFLob0o1Sj1MY61dRUcE3v/lNtm3bxnPPPcf999/Ptm3bctpG79692bBhA5s3b+all15i3bp1PPfcczltQ4qbkliIGhsbWbNmDTfddFOhuyJdKPU4hbV+YRSxDauYsRQvJbEQ3Xbbbdx333306KFhL2alHqeO1i/fxX/zWZy3o2LGUj5K86+0CD366KMMHDiQSZMmFbor0oVSj1Nn67dgwYK8lRrKd3HetmLGjY2NbNq0iS1btuS8DSleSmIh2bhxI6tXr2bEiBHMnj2bDRs2cO211xa6W9JOqccp7PULszhvKRQzlgxkWq+qGCeKuCZfsl/+8pd+6aWXdvo+Eav3lu6kOBXHuHe1frmIUTwe9+uuu85vvfXWvLWxf/9+f+utt9zd/d133/Vzzz3XH3nkkfZtFHzsNeVv0p6YiOTlnNjGjRtZvnw5GzZsoL6+nvr6etauXZvTNvbu3cvUqVOZMGECkydPZtq0aSV7W4R0TDc7F6FSv0FTcSoM3ewspUh7YiIiEllKYiIiElmlVgA4reKZxSqbYphRoDgVRjrjrsKyEhUldU5MRETKiw4niohIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZCmJiYhIZP1/TWJTDs+LG+4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from fealpy.mesh.mesh_tools import unique_row\n",
    "\n",
    "node =np.array(\n",
    "    [(0.0, 0.0),\n",
    "     (1.0, 0.0),\n",
    "     (1.0, 1.0),\n",
    "     (0.0, 1.0)], dtype=np.float)\n",
    "\n",
    "cell = np.array([\n",
    "        (1, 2, 0), \n",
    "        (3, 0, 2)], dtype=np.int)\n",
    "localEdge = np.array(\n",
    "    [(1, 2),\n",
    "     (2, 0),\n",
    "     (0, 1)], dtype=np.int)\n",
    "\n",
    "totalEdge = cell[:, localEdge].reshape(-1, 2)\n",
    "stotalEdge = np.sort(totalEdge, axis=1)\n",
    "uedge, i0, j = unique_row(stotalEdge)\n",
    "NE = i0.shape[0]\n",
    "\n",
    "i1 = np.zeros(NE, dtype=np.int) \n",
    "NC = cell.shape[0]\n",
    "i1[j] = np.arange(3*NC)\n",
    "\n",
    "edge2cell = np.zeros((NE, 4), dtype=np.int)\n",
    "t0 = i0//3\n",
    "t1 = i1//3\n",
    "k0 = i0%3\n",
    "k1 = i1%3\n",
    "edge2cell[:, 0] = t0\n",
    "edge2cell[:, 1] = t1\n",
    "edge2cell[:, 2] = k0\n",
    "edge2cell[:, 3] = k1\n",
    "\n",
    "edge = totalEdge[i0]\n",
    "\n",
    "from fealpy.mesh.TriangleMesh import TriangleMesh\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "tmesh = TriangleMesh(point, cell)\n",
    "fig, axes = plt.subplots(2, 4)\n",
    "tmesh.add_plot(axes[0, 0])\n",
    "tmesh.find_node(axes[0, 0], showindex=True, markersize=25, fontsize=12)\n",
    "tmesh.find_cell(axes[0, 0], showindex=True, markersize=100, fontsize=12)\n",
    "axes[0, 0].set_title('mesh')\n",
    "\n",
    "for ax in axes.reshape(-1)[1:]:\n",
    "    ax.axis('tight')\n",
    "    ax.axis('off')\n",
    "axes[0, 1].table(cellText=cell, rowLabels=['0:', '1:'], loc='center')\n",
    "axes[0, 1].set_title('cell', y=0.6)\n",
    "axes[0, 2].table(cellText=totalEdge, rowLabels=['0:', '1:', '2:', '3:', '4:', '5:'], loc='center')\n",
    "axes[0, 2].set_title('totalEdge', y=0.85)\n",
    "axes[0, 3].table(cellText=np.sort(totalEdge, axis=1), rowLabels=['0:', '1:', '2:', '3:', '4:', '5:'], loc='center')\n",
    "axes[0, 3].set_title('stotalEdge', y=0.85)\n",
    "axes[1, 0].table(cellText=i0.reshape(-1,1),  loc='center')\n",
    "axes[1, 0].set_title('i0', y=0.85)\n",
    "\n",
    "axes[1, 1].table(cellText=i1.reshape(-1,1),  loc='center')\n",
    "axes[1, 1].set_title('i1', y=0.85)\n",
    "\n",
    "axes[1, 2].table(cellText=edge,  rowLabels=['0:', '1:', '2:', '3:', '4:'], loc='center')\n",
    "axes[1, 2].set_title('edge', y=0.85)\n",
    "\n",
    "axes[1, 3].table(cellText=localEdge,  rowLabels=['0:', '1:', '2:'], loc='center')\n",
    "axes[1, 3].set_title('localEdge', y=0.85)\n",
    "\n",
    "#axes[1, 1].table(cellText=i1.reshape(-1,1),  loc='center')\n",
    "#axes[1, 1].set_title('i1', y=0.85)\n",
    "#axes[1, 2].table(cellText=edge2cell,  rowLabels=['0:', '1:', '2:', '3:', '4:'], loc='center')\n",
    "#axes[1, 2].set_title('edge2cell', y=0.85)\n",
    "#axes[1, 3].table(cellText=edge,  rowLabels=['0:', '1:', '2:', '3:', '4:'], loc='center')\n",
    "#axes[1, 3].set_title('edge', y=0.85)\n",
    "plt.tight_layout(pad=1, w_pad=1, h_pad=1.0)\n",
    "plt.savefig('/home/why/edgedata0.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面给出网格的数据结构示意图："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAGFCAYAAADdO6Z6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8VNX9//HXJyFAQKAiEBZFRJTF1iAoGmqRTdGfrVip\ntLUWaN1oQLFg9QsUFOvyRW2pVJaKgvpF60JFbIkiokhUUAmQFtmqgKKSSFxQE5bInN8fd5ImIesk\nM3eW9/PxmAfMnXPnfs7MOckn5557rjnnEBEREZG6SfI7ABEREZFYpCRKREREJARKokRERERCoCRK\nREREJARKokRERERCoCRKREREJARKokRERERCoCRKREREJARKokRERERCoCRKjmJm55lZwMwu8zsW\nERGRaKUkSqqi+wGJiIhUQ0mUVMX8DkBERCSaKYkSERERCYGSqChiZrcF5yKdYmaLzexLM/vUzG4P\nvn6CmT1nZvvNbK+ZTaywf2Mzm2Fm/zGzg2b2oZnNNLPGFcqdb2bZZvaFmX1tZtvM7M4K4Tggycym\nmtkeMztgZi+b2cnh/RRERERiQyO/A5BySuYhPQVsAW4BLgammtnnwHXAKuBm4BfAvWb2tnPudTMz\n4B9Af+CvwDbge8BvgVOAywDMrFew3CZgGnAI6BbcrywDJgNHgHuBVsF4FgMZDV1xERGRWKMkKjqt\nc85lApjZAmA3cB/wP865+4LbnwQ+AX4NvI6XVA0GBjjn1pa8kZm9C8wzs3Occ+uA84EU4CLn3Bc1\nxNEESHfOHQm+15fAn82sl3NuS4PVVkREJAbpdF70ccDDpU+cCwDr8UaGFpbZvh/YDnQNbvoJsBXY\nYWbHlTyAV4P7DgqW+zL474+Do1fVWViSQAVlB9+raxXlRUREEoaSqOj0YYXn+4GDzrnPK9l+bPD/\npwCnAfsqPLbjJWbtguWeAt4AFgD5ZvY3M7u8ioRqT4XnJSNXx1YsKCIikmh0Oi86HanlNvjvUgRJ\nwL/x5kBVmRA55w4CA8xsEN58qwuBnwKrzOwC51zZ9aFqOqZITDOz1YBzzg2qqaxIojOzAHCbc67k\nYqcxeGdIujjnKv7xnxCURMWP94HTnXOv1qZwsNyrwE1mNhm4A++U3yvhC1Ek6jggEMqOZvZzoJ1z\n7v6GDUkkZjgSfGFmnc6LH08Dx5vZNRVfMLOmZtYs+P/KTsXl4o0uNQlviCJR53xgWIj7XgFMaMBY\nRCTGaCQqfvwfMBLvSrxBePOekoGewOXABcAGYLqZDQCWAx8AacBv8OZhve5D3CK+cc5963cMIhK7\nNBIVO6oaMnXgTeoAhgP/A3wXb22n6UBfYBawI1h+GV7y9CvgAbwEajUwxDn3dW2PJ+KHMgvSdjez\np4MLzxaY2Z/NrEmZcslmNs3M3gsuPLvLzO6sZOHZ1Wb2SpnnJTffvry6hWbN7FW8OYUnBssHzGxn\nJD4Dkboys45m9rCZfRzsDzvNbK6ZNQq+3irYhz4Mvv4fM7u5FldwJzyNREUR59wMYEYl23+Fl/RU\n3D6owvMjeOtJ3VfNMVbjJU3VxfEa3ihWxe0fVLZdJIJKkvingV14fzScA9wAfAcYE3z9YWBUsNx9\nwNl4i8f2AEZU8n4V/Q/VLzR7R3B7J+BGvNPh39SnYiLhYGYdgHeAlngLMW/Ha7c/AZqZWTGwBugA\nzMe7CKk/cDfQHphYydtKkJIoEYlF7zvnLgv+f56ZfQ38xszuwxthHwU86JwbGywz38z2AZPM7Lzg\nHwrVqXahWefcKjP7GPiOc+5vDV47kYbzv3hL3PRzzm0ss/02ADP7PXAS0Ns5VzKausDM9uJdePRH\n59zHkQw4luh0nojEGgfMqbDtL3ijQf8v+HB4p7HL+mOwzMW1OIYWmpWYFzwdNxx4vkICVdZP8Nr3\n/goLNa/CG2gZEJloY5NGokQkFr1X4fn7eEsVdMFLdgIVyzjn8oMjSifW4v210KzEg7Z4p/HerabM\nKXj3Wd1XyWtlF2qWSiiJEpF4UNncpvpcBKGFZiVRJAErgZlU3r53VLJNgpREiUgsOgXvKtMS3fB+\nGewO/psULLO9pICZtcObfF52v/rQlaoS7fYBX+FdsV2V94FjartQs5SnOVEiEmsMGFdh2w14Sc0L\nQFawzI0VykwKllneQHEU4l2hJxKVgkvfPAf8yMz6VFHsaSDDzC6o+EJw6QNdkV2NhBqJMrNeeFck\n9MW7dLMI2ALc65z7p4+hRZyZTQX+AGx2zp3udzzhZGZn4l36PhBvzsxnwDrg9865//gWmNTHSWa2\nDHgR73LsXwCLnXP/BjCzR4Frgyv0v4a3xMEo4NlaXJlXWznASDP7I94l5N8k2s8RiQlT8FbmX2Nm\nDwJbgY54E8q/j7eMxyXAP83sEbx23Rw4HbgM72fm5xGPOkYkVBKFN6H0GOAR4BOgGd6aMc+b2bXO\nuYd8jC1izKwT3po5ibKuzS14v2ifAf6Fl0BfD2wws7Odc1v8DE7qzOHdNPsPeGvZfAvMBm4uU+Yq\nvNMUY4BLgTzgTuD2Kt6vuudVbZ8LpAePcSPeaUIlURJVnHOfmNnZeP3lCryJ5h/jjdgWOee+Dd7F\nYgre3S1+iXcKcAfegs37y74dOo1djnmjfYkreAnoBqCJc66X3/FEgpk9CRyHl0QflwAjUecA68ve\n4sPMugH/Bp5xzo3yLTipEzO7Fe8He1vnnP46FhFfJfycqOA54z14E07jXvAvjss4er5I3HLOrat4\njzTn3Ht4l/329CcqERGJdYl2Og8AM2sGpOJNCh0OXATE/arDZpaEd9pjgXPuXd0WiTRgs99BiIhI\nbErIJApv5eLrgv8PAH/HmyMT734DdAYG+x2I38zsSrz7R/3e71hERCQ2JeScKDM7FTge7wqFkcBh\nINM596mvgYWRmbXGmyh4h3Puz8Ftr5IAc6IqMrMeeFfn/RsY4BKxE4iISL0lZBJVkZmtAFo5587x\nO5ZwMbN5eCNQp5XMD0rEJMrM0oA38eYDZjjn8nwOSUREYlTCTywPWgKcZWan+B1IOASvRLsGbz5U\nJzM70cy6AE2BlODzuL8nmJm1xFtXqCVwoRIoERGpDyVRntTgv/G6+nAnvBWcZwO7go+deAsQdg/+\nf5pv0UWAmTXBW8OnG3Cxc257DbuIiIhUK6EmlptZW+fcvgrbGgGjgQN4q5fHo83AjyvZfife4qM3\n4CVScSl4VeLTeEnjJc65t30OSURE4kBCzYkys2fxTuWswVuxtT3e7SK6AxOdc/f7GF7EJcqcKDP7\nM16i+DzequXlOOcej3hQIiIS8xIniTJzTwMPA6uBwwaNmjSiRVpzOp3RgeO6xv6UoMOFh1n71xz4\nOV5aWJVtwJPQsuMxHDkc4MxR6RGKMLy2LN9Brx2fsbpsk3bOgsnigKr2c87pBpsiIlJnCXU6byTw\nI6BDstHvN2cxdMoP/A6pQRUfKGbdQxtwn9eQGH8Blmxcn/1rUlJTIhNcBBw5fITc/3zOAedKJ7kB\nOOcG+RaUiIjErYSaWH4Ab0XNrx2c8fPv+h1Og0tJTaH7sJNJ2mTeEqKVCUDSJqPHhd3iKoECaNej\nDfsDjuvxvmsREZFwSpgkaqB5I1CLkowf/fECjjsp9k/fVSbjmr4EPnWwgqMTqQDwIgQ+dWRc09eH\n6MJn97qPeP6ml0hObcSiJKNjkjEQjvgdl4iIxK+EmRPVrkcb1/38kznj59+N2wSqxDuP5bJ88ssk\ntTMCvR0cC3zhjUAFPnVcfPdQzoqTeVDgJVCLr/g7jZo24vrsX3Hwq0Os+fM6cp/Z8qRz7ud+xyci\nIvEpYZKogZP6u0E39fc7jIj58O2PWbsgh20vvoc74rBk7xRexjV96dyvk9/hNZiKCVTz45oB8PGm\nPB68aHHC32FZRETCJ6EmlieSzv060blfJ4oPFHPom8M0OaZx3M2BqiqBEhERiQQlUXEuJTUl7pIn\nUAIlIiL+S5iJ5RI/lECJiEg0UBIlMUUJlIiIRAslURIzlECJiEg0URIlMUEJlIiIRBslURL1lECJ\niEg0UhIlUU0JlIiIRCslURK1lECJiEg0UxIlUUkJlIiIRDslURJ1lECJiEgsUBIlUUUJlIiIxAol\nURI1lECJiEgsURIlUUEJlIiIxBolUeI7JVAiIhKLlESJr5RAiYhIrFISJb5RAiUiIrFMSZT4QgmU\niIjEOiVREnFKoEREJB4oiZKIUgIlIiLxQkmURIwSKBERiSdKoiQilECJiEi8URIlYacESkRE4lEj\nvwOQyh0uKuaNOW/z8aY8Pt6Ux4EvD3Lpny+k9+Wn+R1anSiBEhGReKUkKkoVfX6A1/68ju8c35L2\np7Vj95t7/A6pzpRAhZ+ZdQba1OMtGgOHGyicWJPIdS9wzn3odxAisU5JVJRqkdacm3J/wzFtmvHJ\nv/J58KLFfodUJ0qgws/MOiclJW0PBAJNQ32PpKQkAoFAQ4YVMxK87gfNrLsSKZH6URIVpZJTkjmm\nTWwmHkqgIqZNIBBounjxYnr27FnnnbOyspg2bRqh7h/LErnuW7du5corr2yKN4KpJEqkHpRESYNS\nAhV5PXv2pE+fPnXeb+vWrfXaP5Ylct1FpOHo6rwo4Zyj8LMivtizn8LPinDO+R1SnSmBig5z5szh\npJNOIjU1lXPOOYd33nnH75DqrK51eOaZZ+jZsyepqamkp6fzwgsvRCjShpfIdReJNUqifHZg/0HW\nPbSBOQMe5t7vzeP+sx/i3u/NY86Ah1n30AYO7D/od4i1ogQqOjz11FNMmjSJGTNmsHHjRtLT0xk2\nbBgFBQV+h1Zrda3Dm2++yRVXXME111zDpk2bGD58OJdeeilbtmyJcOT1l8h1F4lFSqJ89N7q3dx/\n1oOsnPEq3+28n1tvhfvug1tvhe923s/KGa9y/1kPsmf9J36HWi0lUNFj1qxZXHfddYwaNYoePXow\nf/58mjVrxsKFC/0OrdbqWofZs2dz0UUXMXHiRLp3787tt99Onz59eOCBByIcef0lct1FYpGSKJ+8\nt3o3T4x6lvTvFvP003DrdBg4EPr29f69dTo8/TSkf7eYF299xe9wq6QEKnoUFxeTk5PDkCFDSreZ\nGUOHDmXt2rUAjBkzhkGDBvkVYo1qU4eK1q5dy9ChQ8ttGzZsWJXlo1Ui110kVimJ8sGB/QdZcu3z\nnHWm4847oHXrysu1bg133gGn9QQcFBcVRzTOmiiBii4FBQUcOXKEtLS0ctvT0tLIy8sDoGPHjpx4\n4ol+hFcrtalDRXl5eXUqH60Sue4isUpJlA9yn9lC8YFibv4dJCdXXzY5GX75S+//H74TPaf1lEDF\nprvuuotHHnnE7zBEROKCljiIMOcc6x/dwIABVY9AlVi6FAoLYd8+7/n2Ff/htT+vwwzOvqoPTY5p\nHP6AK6EEKjq1adOG5ORk8vPzy23Pz8+nffv2PkVVN6HUoX379jFd5xKJXHeRWKWRqAgr+vwABe/v\nZ8CAmss+/TQsWgT//CeYweHCb1l935u8eu+bHPjSn6v2lEBFr5SUFPr27cuqVatKtznnWLVqFf37\n9/cxstoLpQ4ZGRnlygOsXLmSjIyMsMba0BK57iKxSiNREXY4OK+pRYuay/7tb//9//r18LvfwQ1r\nr+LYE1qFKbrqKYGKfhMnTmTMmDH07duXfv36MWvWLIqKihgzZgwAkydP5pNPPuHRRx/1N9Bq1FSH\nUaNGcfzxx3PXXXcBMGHCBAYOHMif/vQnLr74Yv72t7+Rk5PDggULfKxFaBK57iKxSElUhDVulgLA\n11/Xbb+S8k2a6xSeVG3kyJEUFBQwffp08vPz6d27NytWrKBt27aANxF5z57ovpl1TXX46KOPaNTo\nvz+6MjIyeOKJJ5g6dSpTp07llFNOYdmyZfTq1cuvKoQskesuEouUREVYs9aptDm5Fa+t2c/AgbXf\nb80aaHNyK1KPDflesyFTAhVbMjMzyczMrPS1RYsWRTia0FRXh1deOXrJjxEjRjBixIhwhxURiVx3\nkVijOVERZmacOboP2Wvg889rt89nn0F2Npw5ug9mFt4AK1ACJSIiUjmNRPkg/fJerL7nde65t5g7\n76h+mYMjR+De+4yU1EakXx7ZIXolULEhKyur9Ia6dfHGG2/Ua/9Ylsh137Vrl98hRISZdQba1OMt\nGgOHGyiccFKcDavAOfdhbQtbLN7oNhQDJ/V3g26KniuUSlYsP+tMx+9uguOOO7rMZ5/BvffBO+uN\nKx67jG4Du0QsvlhPoD7elMeDFy2O7LBdhJnZOUlJSWsDgUDI75GUlER99o9lqnsgwzm3zu9YwsHM\nOiclJW0PBAIhz3+IlfahOBtWUlLSwUAg0L22iZRGonzSbWAXrnjsMpZc+zw//WkxP/gBDBjgXbX3\n9dfeHKjsbEhJTeGKxy5RAiWVORwIBFi8eDE9e/as885ZWVlMmzYt5P1jWSLXfevWrVx55ZUQG6MC\noWoTCASaxnvfUJwNK9g3muKNYCqJinbdBnZhwjvXkvvMFtY/uoHVt+8vfa3Nya04/9Y+9B55Gk1b\nNolYTEqgYk/Pnj3p06cP2dnZ3HvvveTk5LB3716ee+45Lrnkkir3KzmNVbJ/NKhrHQBWr17NpEmT\nePfdd+ncuTNTp05l9OjR1e4TjXUvMWfOHO677z7y8vJIT0/nL3/5C2eddVaV5Z955hmmT5/O7t27\nOfXUU/nf//1fLrrooghGHL3Kfr91+Vz9bB916QN+xnn33XezdOlStm3bRmpqKv3792fmzJmceuqp\nURUn1L1P1UVCTSz/ODeP5VNWMWfQI9zZbTazznqQZ677B5/t/MK3mFJbNeWcq/swbs1V3Lw5kwlv\nXc3NmzMZt+Yqzrm6T4MmUIeLinn13jdY/Iu/M/O0OdzW6Y9seubd0tfjNYEys15m9rSZvW9mhWa2\nz8xeM7Mf+h1bQyosLKR3797MnTs34hcgNJS61mH37t388Ic/ZMiQIeTm5jJhwgSuvvpqVq5cGYFo\nG95TTz3FpEmTmDFjBhs3biQ9PZ1hw4ZRUFBQafk333yTK664gmuuuYZNmzYxfPhwLr30UrZs2RLh\nyKNbXT9XP8VKP87Ozub666/nrbfe4uWXX6a4uJgLLriAAwcO+B1aOeH+7hNqJOqNOe+wZ/0n9Prh\nqaT1bMM3+4p4e+FG/jrs/7h6+S9od2olE5MixMxo1jqVZq1Tw3aMos8P8Nqf1/Gd41vS/rR27H7z\nv+sFxWsCFXQicAzwCPAJ0AwYATxvZtc65x7yMbYGc+GFF3LhhRcC3krXsaiudZg3bx5du3blnnvu\nAaB79+68/vrrzJo1i/PPPz+ssYbDrFmzuO666xg1ahQA8+fPZ/ny5SxcuJCbb775qPKzZ8/moosu\nYuLEiQDcfvvtrFy5kgceeIC5c+dGNPZoVtfP1U+x0o+zsrLKPX/kkUdo164dOTk5nHvuuT5FdbRw\nf/cJNRKVcV1fbnz7Gi66fRB9fv49BtxwNr9a+lMCRxyvP/C23+GFXYu05tyU+xtufOsazv/9gNIO\nGucJFM65F5xz/8859wfn3MPOub8Ag4BcYKLP4UXMbbfdxkknneR3GA1q3bp1DB06tNy2YcOGsXbt\nWp8iCl1xcTE5OTkMGTKkdJuZMXTo0Crrs3bt2ripf7jU5nMdM2YMgwYN8ivEuPDll19iZrSu6aaw\nERRKn6qrhEqiTujbkeRG5at83EnH0vbU4yj4z2c+RRU5ySnJHNOmfIJU8N7ncZ1AVcV5GeQe4Dt+\nxxIpbdu2pVu3bn6H0aDy8vJIS0srty0tLY2vvvqKQ4cO+RRVaAoKCjhy5Eil9cnLy6t0n6rqX1X5\nRFSbz7Vjx46ceOKJfoQXF5xz3HjjjZx77rlRtVp+KH2qrhLqdF5VCvcV0q5HfZYTiV1vzl9P4+aN\nEyKBMrNmQCrQChgOXAT8rdqd4si4ceMYN26c32GIRJ2SexFKaDIzM9myZUvp+muJJKFGoiqT+/ct\nfJX3Dd8d3sPvUCLqk3/ng/NGpxIhgQr6I7APeA+4F3gWuN7XiKRe2rdvT35+frlt+fn5tGzZkiZN\nIndVa0No06YNycnJldanffv2le5TVf2rKp+IQvlcpfbGjx9PVlYWq1evpkOHDn6HU04kvvuETqL2\n/eczsqa+QuezOkV8NXA/7V73ES9M9e7BNXTqDxIlgQKYBQwFRgFZQDIQW79ppZyMjAxWrVpVbttL\nL71ERkaGTxGFLiUlhb59+5arj3OOVatW0b9/5QsFV1b/lStXxmT9wyWUz1VqZ/z48SxbtoxXX32V\nzp07+x3OUSLx3SdMEnW48DDFB4pLn3+zr5AnRi0ltVUTLn/wR1F9KWl9FB8o5pt9haV1L5lEntw4\nGQyaHNPY5wjD49tD32Jm5U6EO+d2OOdecc4tds5dArQAnvcnwoZXWFhIbm4umzZtAmDnzp3k5uay\nZ493FeacOXOOmoQcbWqqw+TJk8utATV27Fh27tzJLbfcwvbt25k7dy5LliwpvVot1kycOJEFCxbw\n2GOPsW3bNsaOHUtRURFjxowBYNSoUUyZMqW0/IQJE3jxxRf505/+xPbt27ntttvIyclh/PjxPtUg\nOtX0uVZsV36qqQ9Ei8zMTB5//HGeeOIJmjdvTn5+Pvn5+Rw8eNDv0Mqp6buvr4SZE7X2rzmse2gD\n3YedzJlXpvPynWs49PVhfv3cz2jRrrnf4TW4D976iLULcti+4n3cEYclGyec1ZGPN+SR0iyFy+f/\nkP/7+RK/w2xwpfV+8T2APKC67HgJMN/MTnHO/SciAYbR+vXrGTRoEGaGmTFp0iQARo8ezcKFCyko\nKGDnzp0+R1m9muqQl5dX7pdJly5dWL58Ob/97W+ZPXs2xx9/PA8//HDUJ4tVGTlyJAUFBUyfPp38\n/Hx69+7NihUraNu2LQAfffQRjRr998d2RkYGTzzxBFOnTmXq1KmccsopLFu2LKom90aDmj7Xiu3K\nTzX1gWgxf/58zIyBAweW275o0aLS5QSiQU3ffX0lTBLFz8F97ti+/j22Zb1HcuNkfrX0p7TpFj2X\nYzaUdx7dxPIpq0hqZ7ihDlp7df9w/cdwGAbf8n1SvxPyLaWiVrl6nw/U/NWWLMrVKqyBRch5551X\n7b2pbr31Vm699dYIRlR3NdVh0aJFR20bMGAAOTk54QwrojIzM8nMzKz0tVdeeeWobSNGjGDEiBHh\nDivmVfe5Vtau/FJTH4gWsRBjieq++/pKnCSqOxAAtwvYB0cOHyFQHDuNoLY+eOsjlk9ZBf0gMMyV\nP2F7NvAirLxjDZfce4FfIYZFdfU2s7bOuX1ly5tZI2A0cADQ8s4iIlJniZNEAawAdgDdwfbCi7e9\nyjlX9y19+fTLovfGiLW1dkEOSe3s6AQKvOfHgTWDNbPfAmD7S+/z1SdfA3D2VX1ido5UtfWGv5pZ\nS2AN8DHQHvgFXmo90TlXFNFgRUQkLiRWEpWPN0NmBzgHn2zKZ+kNL5S+HOtJVPGBYm8O1NBKEwnP\nWnCF8GXhfizJ2PbCe2x74T0ATh/RKyaTqFrU+0ngKmAscBzwNZAD/M45tzxigYZJVlZW6Q0+66Jk\nTZdQ949liVz3Xbt2+R1CxMR731CcDSuUvmHRfG+ehmQzrHxFtwFPwk25YzmmbXxMLP9mXyH3pc+H\nn+ONsVQlzupeXb3drS4+L7sEzOycpKSktfWZm5CUlBRTcxsakuoeyHDOrfM7lnBIpL6hOBtWXftG\nYo1ElfUFWLLF5MhLVZoc0xhLNtznNSTGcVb3Wtc7/hwOBAIsXryYnj3rPoqalZXFtGnTQt4/liVy\n3bdu3cqVV14JcNjvWMIoIfqG4mxYofSNxEyiApC0yeh+YTdSUlP8jqbBpKSmcMJZHb2r8M6m8lNb\ncVj3lNQUug87mR057xM4u5pTmXGqZ8+e9OnTh7vvvpulS5eybds2UlNT6d+/PzNnzuTUU0+tdL+S\nYfWS/aNBdnY29957Lzk5Oezdu5fnnnuOSy65pNp9Vq9ezaRJk3j33Xfp3LkzU6dOrXHNn0SueyIp\n+X7r+tn62T7mz5/PvHnz2L17NwCnnXYa06dP58ILLzyqbDS04zlz5nDfffeRl5dHeno6f/nLXzjr\nrLPKlfEzzrr+XKyrBPt1AwSAFyHwqSPjmr41Fo8lu9d9xMcb8uAzvEn0FUdO47juGdf0JfCpq7ze\nCSI7O5vrr7+et956i5dffpni4mIuuOACDhw44HdotVZYWEjv3r2ZO3durRbA3b17Nz/84Q8ZMmQI\nubm5TJgwgauvvpqVK1dGINqGlch1D7e6frZ+OuGEE5g5cyYbNmwgJyeHwYMHM3z48KicS/TUU08x\nadIkZsyYwcaNG0lPT2fYsGEUFBT4HVqpcP9cTJyRqG3AF94oTOBTx8V3D6Vzv05+R9VgSlYiT2mW\nwuBbvs/KO9aQtNsI9HZwLHFdd4ATzz6ei+8eyvLJL5evdwLJysoq9/yRRx6hXbt25OTkcO655/oU\nVd1ceOGFpX9x12a+5rx58+jatSv33HMPAN27d+f1119n1qxZnH/++WGNtaElct3Dra6frZ8uvvji\ncs/vuOMO5s2bx7p166LuVNisWbO47rrrShfXnD9/PsuXL2fhwoXcfPPNPkfnCffPxcQZiXoS7GWj\n+5nd+PXSn3HWqHS/I2owJQlUo6aNuD77V3z/N2fx66U/o/uZ3bCXLa7rXtZZo9L/W++VeNfkJbAv\nv/wSM6N16/+uOvqrX/2KQYMG+RhVw1q3bt1Rq5MPGzaMtWvX+hRR5CRy3RvabbfdxkknneR3GEcJ\nBAI8+eSTFBUVRd39EIuLi8nJyWHIkCGl28yMoUOHRnUbrOznYn0kzEhUxnV9GXzz9+NmHlCJiglU\nyc2EO/frROd+nSg+UMyhbw7T5JjGcVf3ypTU+4O3PmLRj59K2Fu0O+e48cYbOffcc8vdAqRDhw5R\n/5d4XeTl5ZGWVu4WiaSlpfHVV19x6NAhmjSJ3/tLJ3LdG1rbtm3p1q2b32GU2rx5MxkZGRw8eJAW\nLVqwdOkgFssiAAAgAElEQVRSevTo4XdY5RQUFHDkyJFK2+D27dt9iqp6Vf1crI+ESaIaN4+/JKKq\nBKqslNSUuKt3bTRq0gjnXL7fcfglMzOTLVu2lK7PUuKuu+7yKSKR6DVu3DjGjRvndxilevToQW5u\nLvv372fJkiWMGjWKNWvWRF0iFWuq+rlYHwmTRMWb2iRQkfDt4SO8es8b/OvZrRz48iBpvdoy+Obv\nc/KAE32JR2D8+PFkZWWRnZ1Nhw4d/A4nrNq3b09+fvlcOT8/n5YtW8b9SEwi1z3eNWrUiK5duwJw\nxhln8Pbbb3P//fczb948nyP7rzZt2pCcnFxpG2zfPvpOAoTr52LizImKI9GSQAE8N+EF1j20gdNH\n9OSiPwwiKdl4/JfP8uE7H/sWUyIbP348y5Yt49VXX6Vz585+hxN2GRkZrFq1qty2l156Kermj4RD\nItc90QQCAQ4dOuR3GOWkpKTQt2/fcm3QOceqVavo37+/j5EdLZw/F5VExZhoSqA+2riXzc9vZ+iU\nH3D+1AH0/cXpjH76cr5zfEtW3rHGt7gSVWZmJo8//jhPPPEEzZs3Jz8/n/z8fA4ePFhaZsqUKVG9\njlBhYSG5ubls2rQJgJ07d5Kbm8uePXsAmDx5crn4x44dy86dO7nlllvYvn07c+fOZcmSJUycONGX\n+OsjkesebjV9tnPmzDlqkr5fpkyZQnZ2Nh988AGbN29m8uTJvPbaayWLQEaViRMnsmDBAh577DG2\nbdvG2LFjKSoqYsyYMX6HVqo2PxfrQ6fzYkg0JVAAW/65g6RGSfT5xfdKtzVq0ogzfvZdXpn5Bl/t\n/ZqWHVr4GGFimT9/PmbGwIEDy21ftGhR6SXIe/fuLf3FEY3Wr1/PoEGDMDPMjEmTJgEwevRoFi5c\nSF5eXrn4u3TpwvLly/ntb3/L7NmzOf7443n44Yej5hdiXSRy3cOtps+2oKCAnTt3+hyl59NPP2X0\n6NHs3buXVq1acfrpp/PSSy8xePBgv0M7ysiRIykoKGD69Onk5+fTu3dvVqxYQdu2bf0OrVRtfi7W\nh5KoGBFtCRRA3rv7OK7rsTRpXv72MZ3O6FD6upKoyKnNfakWLVoUgUhCd95551Vbj8riHzBgADk5\nOeEMKyISue7hVtNne+utt3LrrbdGMKKqPfTQQ36HUCeZmZlkZmb6HUaVwn2/PiVRMSCaEijnHEWf\nH+BwUTFfffI1LTscc1SZFmnNcc7xdf43PkQoIiISGUqioly0JFAH9h8k95ktrH90AwXv7y/dvv+T\nr1j30AbSL+9FaqumgHdKD6D4wLe+xJposrKyQrolRMllvqHuH8sSue67du3yO4SICfW7LfmMor1t\nKM6GFUp8Fk8L71Vn4KT+btBN0XXFQE2iJYF6b/Vullz7PMUHivnBADhvALRoAffcA8XF8NVX3npU\nP3nwEroN7MK+/3zGnIGP8KN7zqfvL073JeaPN+Xx4EWLo/smWfVkZuckJSWtrc9wdVJSUtiHu6OV\n6h7IcM6t8zuWcDCzzklJSdsDgUDTUN8jVtqH4mxYSUlJBwOBQHfn3Ie1Ka+RqCgVTQnUE6Oe5awz\nHTf/DsqulN+5MxQUwEMPwT33FvPEqGe54rHLSGrkXfTZIu3oU33SoA4HAgEWL14c0j21srKymDZt\nWsj7x7JErvvWrVtLrvQ67Hcs4eKc+9DMugNtQn2PQCDQmBj4jBRnwwoEAgW1TaBASVRUipYE6sD+\ngyy59nnOOtNx5x2QnFz+9ZNPhk2bIDUV7rwDpv7eseTa5+l3dR/MjPanRc8VGvGsZ8+e9OnTp877\nlQxdh7p/LEvkuieK4C/CWv8yFAmF1omKMtGSQAHkPrOF4gPF3Py7oxMogPPOgyNH4B//8F7/3U1w\nuKiY9f/3L47v00FX5kXY/PnzSU9Pp1WrVrRq1Yr+/fvz4osv+h1WSObMmcNJJ51Eamoq55xzDu+8\n80615Z955hl69uxJamoq6enpvPDCCxGKtOFkZ2dzySWX0KlTJ5KSknj++edr3Gf16tX07duXpk2b\ncuqpp/Loo49GIFIRKaEkKopEUwLlnGP9oxsYMKD8Kbyyevb0EqmHHoK//hXWroVjjoGizw4w9Pc/\niGzAwgknnMDMmTPZsGEDOTk5DB48mOHDh0f9ZM6KnnrqKSZNmsSMGTPYuHEj6enpDBs2jIKCgkrL\nv/nmm1xxxRVcc801bNq0ieHDh3PppZeyZcuWCEdeP4WFhfTu3Zu5c+diVvN0vt27d/PDH/6QIUOG\nkJuby4QJE7j66qtZuXJlBKIVEdDpvKgRTQkUQNHnByh4fz8Dalgkd8oUWLgQXn4Zvv4a2rTxJpq3\n7XZcZAKVUhdffHG553fccQfz5s1j3bp1MTXvZ9asWVx33XWlC+HNnz+f5cuXs3DhQm6++eajys+e\nPZuLLrqodKXu22+/nZUrV/LAAw8wd+7ciMZeHxdeeCEXXngh4P0RU5N58+bRtWtX7rnnHgC6d+/O\n66+/zqxZszj//PPDGquIeDQSFQWiLYEC77QceFfhVSclBa67Dp55Bl58EW680dt+qDDq5w/GtUAg\nwJNPPklRUVHpfazGjBnDoEGDfI6sesXFxeTk5DBkyJDSbWbG0KFDWbt2baX7rF279qhVuocNG1Zl\n+Xixbt26hKy3SDTRSJTPojGBAmjcLAXwRpfqoqR8xVXMJTI2b95MRkYGBw8epEWLFixdupTu3bsD\n0LFjR5+jq1lBQQFHjhwhLS2t3Pa0tDS2b99e6T55eXmVls/LywtbnNGgqnp/9dVXHDp0iCZNmvgU\nmUjiUBLlo2hNoACatU6lzcmteG3Nfirccqhaa9ZAm5NbkXpsyMuzSD306NGD3Nxc9u/fz5IlSxg1\nahRr1qyhR48e3HXXXX6HJyISV3Q6zyfRnECBdwrlzNF9yF4Dn39eu30++wyys+HM0X1qNTFWGl6j\nRo3o2rUrZ5xxBnfeeSfp6encf//9fodVa23atCE5OZn8/Pxy2/Pz82nfvn2l+7Rv375O5eNFVfVu\n2bKlRqFEIkRJlA+iPYEqkX55L1JSU7jnXm8pg+ocOQL33mekpKaQfnmvyAQoNQoEAhw6dMjvMGot\nJSWFvn37smrVqtJtzjlWrVpVOrerooyMjHLlAVauXElGRkZYY/VbZfV+6aWX4r7eItFESVSExUoC\nBZDaqik/efAS3llvTP29N9JUmc8+g6m/h3fWw08evKT0HnoSWVOmTCE7O5sPPviAzZs3M3nyZF57\n7bWS1amZPHkyo0eP9jnKmk2cOJEFCxbw2GOPsW3bNsaOHUtRURFjxowBYNSoUUyZMqW0/IQJE3jx\nxRf505/+xPbt27ntttvIyclh/PjxPtUgNIWFheTm5rJp0yYAdu7cSW5uLnv27AGO/v7Gjh3Lzp07\nueWWW9i+fTtz585lyZIlpVcpikj4aU5UBMVSAlWi28AuXPHYZSy59nl++tNifvADGBC8d97XX3tz\noLKzvXvnXfGYd+888cenn37K6NGj2bt3L61ateL000/npZdeYvDgwYA3EbnkF3I0GzlyJAUFBUyf\nPp38/Hx69+7NihUraNvWWwH/o48+olGj//7oysjI4IknnmDq1KlMnTqVU045hWXLltGrV2yNiK5f\nv55BgwZhZpgZkyZNAmD06NEsXLjwqO+vS5cuLF++nN/+9rfMnj2b448/nocffvioK/ZEJHyUREVI\nLCZQJboN7MKEd64l95ktrH90A6tv31/6WpuTW3H+rX3oPfI0mrbUPAw/PfTQQ9W+vmjRoghFUn+Z\nmZlkZmZW+torr7xy1LYRI0YwYsSIcIcVVuedd161N2it7PsbMGAAOTk54QxLRKqhJCoCYjmBKpHa\nqinnXN2Hs686gwNfHORQ4WGaNG9M6rFNNYlcREQSkpKoMIuHBKosM6NZ61SatU71OxQJCvW2Lrt2\n7arX/rEskeueiHUWCRclUWEUbwmURJ2CpKSkg1deeWXIM/mTkpJKJ54nmgSv+8FAIFD5zQhFpNaU\nRIWJEigJN+fch2bWHWgT6nsEAoHGQELeoyfB617gnPvQ7zhEYp2SqDBQAiWREvxFqF+GIiI+0DpR\nDUwJlIiISGJQEtWAlECJiIgkDiVRDUQJlIiISGJREtUAlECJiIgkHiVR9aQESkREJDEpiaoHJVAi\nIiKJS0lUiJRAiYiIJDYlUSFQAiUiIiJKoupICZSIiIiAkqg6UQIlIiIiJZRE1ZISKBERESlLSVQt\nKIESERGRipRE1UAJlIiIiFRGSVQ1lECJiIhIVZREVUEJlIiIiFRHSVQllECJiIhITRr5HUCkfXv4\nCK/e8wb/enYrB748SFqvtgy++fucPOBEIL4TqJrqHq/MrDlwM9Av+DgWGOOce8zXwEREJKYl3EjU\ncxNeYN1DGzh9RE8u+sMgkpKNx3/5LB++83FcJ1BQfd3jXBtgGtAD2AQ4f8MREZF4kFAjUR9t3Mvm\n57cz7NaBZFzbF4D0n/Ri7uBH+ectL/P57i/jNoGqru4r71jDVct+7nOEYfUJ0N4596mZ9QXe8Tsg\nERGJfQk1ErXlnztIapREn198r3RboyaN6PL9zny6rYDkxslxmUBB1XU/42ff5aOcvXy192sfowsv\n51yxc+5Tv+MQEZH4klBJVN67+ziu67E0ad64dNvudR+x6anNAFx895C4TKCg8roDdDqjQ+nrIiIi\nUnsJlUR982khLdo1L31edg4UBsVFxT5GF14V616iRVpznHN8nf+ND1GJiIjEroRKoooPfktyE28a\nWNkEavTTl3uvH/jWz/DCovhAMd/sK6S4qLi07mU1Cm6Lx7qLiIiEU0JNLE9p2ogjh7496iq8os8P\neK+nxs/H8cFbH7F2QQ7bV7yPO+JdjHao8DAfvv0xnft1Ki337SEveYqnuouIiERCwoxEbVm+g8OF\nhyl47/OjljH4Or8QgBZpx/gcZcN459FNLLrsKXbkvI8b6uDnwHFQ9OUBFv74Sd55LLe0bLzVHeCz\nnV/w9qKNmNnf/I5FRETiV8IMP/Ta8RnrHBwCkpuWvwrvow17MTPan9bW3yAbwAdvfcTyKaugHwSG\nuf+myR8A64A+sHzyy6T1aEPnfp3iqu4AG5/czD9ueomWwHkwEi+FFBERaXAJMxK12sFLwf8HDh3h\n4FeHAG8V701Pv8vxfTrQskML/wJsIGsX5JDUzmAY5b/dXkAAaA1J7Yy1C3Liru6f7fyCf9z0ElcF\nHJ8EHKsTqH2LiEjkmXMJsnizmQP4CfB3oG334+gxrBs7Vr7Pvv98zsV3D6H9ae38jbGevj30LY+M\neAp3PpBRSYFngG3ACcAH0K77cRS890Vc1B3g7UUb+eDvW9kbcDQt2eicAZjZOOA7QCdgLPAssDFY\narZzLn4XyhIRkbBIuCTqMNAN3B44GHzlX8DvnXMv+xZbAzGzNCCPnwPdKynwLfAqsAHw5tJvAG6J\nh7oDmNnfzoORq8uOQP03idoFdK5i15Occx+GP0IREYkniZNEiYiIiDQgzRkRERERCYGSKBEREZEQ\nKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBER\nEZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGS\nKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERER\nCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKIkS\nERERCYGSKBEREZEQKIkSERERCYGSKBEREZEQKImKQ2Z2npkFzGyA37GI1Ee8t2Uz221mC8s8j+v6\nSngE28z0Ms/HBLd1joZ44pmSqPjl/A5ApIGEtS2bWaqZjTOzFWb2iZl9ZWYbzGysmYX7Z2RldVPf\nlfpylGlH5hljZsvM7EMz+8bM/m1mU82siY9xxjwlUSKS6LoCs4P//yMwCdgJzAUe9isokQbUDFgI\ntAHmAROAt4AZQJaPccW8Rn4HICLiszzgu865rWW2LTCzh4ExZvYH59xOn2ITaQiHgf7OuXVltj1s\nZh8At5nZYOfcKz7FFtM0EhXFzKyjmS00szwzO2hmm83sVxXKdDKz54LDs/lm9iegCWCVvN84M3vf\nzIrMbJ2ZnWtmq83slQrlGpvZDDP7T/C4H5rZTDNrHN4aS7yK5rbsnPusQgJVYmnw354V3rOVmc0y\ns13B99xjZo+aWeu6HFck2C8eNrOPg+1kp5nNNbNGwddbmdmfg+3nYLA93WxmR/WJ6jjniiskUCWW\n4vWvim28Vu03WG6WmX0aPA3+nJl1qqKuA81svZkdCL7vtWZ2m5kFKil7ZbBskZl9ZmZ/M7Pj61Ln\nSNFIVJQys3Z4w61H8E41FAAX4f310MI5N9vMmgKvAMcD9wN7gV8Cg6kwr8LMfgP8BXgN+BPQBXgO\n+ALYU6acAf8A+gN/BbYB3wN+C5wCXBaWCkvciuG23CH4b0GZ92wOvA50xzvVtxHvFMklwdg/Vx+S\n2jCzDsA7QEu8drId6AT8BGhmZsXAGrx2OB+vbfcH7gbaAxMbIIzK2nhd2u/DwBXA48BavP66nKP7\n7BnAC8AnwDS83GNa8LgVy04FbgeeBBYAbYEbgNfM7Azn3Ff1rHPDcs7pEYUP4CHgI+A7FbY/AXyO\n9xf6BLxfTJeVeb0psCO4fUBwWwqwD6+RJ5Up+0sgALxSZtuVQDGQUeG41wbf8xy/Pxs9YusRi205\neJx3gf9UOM6M4L6XVLNvrY8L7AIWlnl+Xtn66hG/D+DRYDs5o4rXfw98BXStsP0uvNNzncpsCwDT\nyzwfHWxHnWuIYSXeHx8ty2yrVfsFTg8ed3aFcouD5crG8zzwNZBWZlvXYD2OlNnWOXjsWyq8Z69g\n2f/x+3ur+NDpvOh1Gd5fA8lmdlzJA3gJaAX0wftrfq9z7tmSnZxzB4EHK7zXmcBxwALnXNmh0yfw\nOlBZPwG2AjsqHPdVvGHfQQ1WQ0kUsdiW5wA9gPEVjnMZkOuce76afdWHpFrB0Z7hwPPOuY1VFPsJ\nkA3sr9COVuGN5NRrGQwzm4I3cnSLKz+6U9v2ezHeKNJfKrz1nylzCt68K1yHAM855/JLtjtvnuEL\nFfYdEdz3mQrH/hTvD5qo6zs6nReFzKwt8B28zP+6Soo4oB1wIvBeJa9vr/D8xOA+75d7E+eOmNnu\nCmVPwfvlsa+a44rUSiy2ZTP7HXA1MNU5t6LCyycDSyrbr77HlYTSFu803rvVlDkF7zRag7cjM/sp\n8AfgIedcxT9Uatt+O+ONRL1foUzFPtsOSKXy/l1xWze8udqVlXV4o1FRRUlUdCoZIVyMN+RbmX+H\n8dj/xjv/XdnkxT2VbBOpSky1ZTMbA/wvMNc5d3ekjitSiSS8020zqbwd7QjlTc3sfLy++A/gN1Uc\n16/2m4SXmF0Y/Leib8J47JAoiYpO+/DOHye7ai47Ne/y1NMqealHhecf4HWGbniTcUv2T8ablJtb\npuz7wOnOuVdDilykvJhpy2Y2HG8i6xLn3Pgqir0PfLeGt1Ifkprsw5vvVF1beh84piHbkZmdDTwL\nvA38tMKp6rLHrU37/QAv6TkZ71RbiYp99lPgIF6freiUSo5twG7nXGWjUVFHc6KiULBh/x0YYWZH\n/WIxszbB/2YBHc1sRJnXmgHXVNhlPfAZcI2VX4H5SuDYCmWfBo43s4rvgZk1Db6/SK3ESls27zYr\nfwNWB9+rKn8H0oMJV1XUh6Razpst/RzwIzPrU0Wxp4EMM7ug4gvBpQ+S63JMM+sJ/BNvIdkfOecO\nVXPc2rTfF/ASnhsqFLuRMlfcBX8GvAxcambty7xXN7wRp7KexRuBurWKOrSubLufLDjzXaJM8LLw\ndXjnkxcAW4DWQF9gsHOujZml4v3lXfGy8EZ4V04Mcs6tCb7fOLzLy1/H6yRdgDEELwt3zg0Jliu5\nvPXCYLk3gGS8dUQuBy5wzm0Ib+0lnkR7Wzbv/mL/Ch7rd3gjBGX9yzn37+B7NsdbruFUYBGQgzfR\n/UfAdc65f9elD5nZLuBV59yvg8/Pw5vAO7CkvhKfzKwj3hIHrfAuoNgKdMSb2P19vKvUsvHa/yN4\nba158PllQBfn3OfB9woAtznnbg8+H423QvlJzrkPzewYvH7XAZiCt9RAWe+74DpSdWy/jwM/w7uw\n4028CeQnA+kV4ukTfP0TvBXTGwHjgHygt3OuNCE0s1vwrkBci5dofo13Jd+lwF+dc3+q62cdVn5f\nHqhH1Q+89WdmA7vxhkM/xrui6ddlyhyPt2Da13gN8o/A+VRymTReo90JFOE10HPwOvHyCuWSgZvw\nfrEU4a3l8TYwFW942ffPRo/YekRzW+a/ywpU9Zhe4T2/g5fofQgcwDut8TBwbF2OGyy3E3i4zHMt\ncZBAj2CbX4S3an4R3mmx+4FGwdebAXfgTdY+EOwX2XijPcll3ucIMK3M83JLHOBdkFFdG19YIa7a\ntt/GwCy8U3ZfBftvx4rxBMsOxBtJPoA3n+tXwL1AYSWfy6V4p+u/Cj7eDX4u3fz+zio+NBKVwIJ/\ncewD/u6cq+zKKZGYoLYsEnvMbCnQyznX3e9YQqU5UQnCKr9T92i80yqaACsxQ21ZJPYE70pQ9vkp\nwP8jxvusRqISRHCuxSzgGbyJuX2BX+MNk57pnPvWx/BEak1tWST2mNkneHO7duLNYxyLd2eAPs65\nimtNxQwtcZA4duPN4bge7y/2z/Ea9GT90pEYsxu1ZZFY8wLeJPT2wCG8ieZTYjmBAo1EiYiIiIRE\nc6JEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRE\nRCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAk\nSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQESqJERERE\nQqAkSkRERCQESqJEREREQqAkSkRERCQESqJEREREQqAkSkRERCQEjfwOoLbMrDPQph5v0Rg43EDh\nxJpErXuBc+5Dv4MIpzjoF34eP5HrHvd9QyQSYiKJMrPOSUlJ2wOBQNNQ3yMpKYlAINCQYcWMRK17\nUlLSQTPrHq+/LOKhX/h5/ASve1z3DZFIiYkkCmgTCASaLl68mJ49e9Z556ysLKZNm0ao+8eyRK37\n1q1bufLKK5vijdLE6y+KmO4Xfh4/keueIH1DJDKcc1H/APoALicnx5V44IEHXJcuXVzTpk3d2Wef\n7d5++21XlcWLF7uK+8equtTbufip+5o1a9yPfvQj17FjR2dmbtmyZdWWz8nJcYAD+rgoaMPheFTW\nL5yrfRsJd9uoKY5wHt+vY9e2nfp5/EToG3roEalHTE4sf+qpp5g0aRIzZsxg48aNpKenM2zYMAoK\nCvwOLawStd4AhYWF9O7dm7lz52JmfocTtaKljfgZh5/H9rud+n18kYTjdxZXmwcV/uI+++yz3Q03\n3OBKBAIB16lTJzdz5kxXmXgZjalrvZ2Ln7qXpZGoyvuFc3VrI+FsG7WJI1zH9/PYZfkxElWb4ydC\n39BDj0g9Ym4kqri4mJycHIYMGVK6zcwYOnQoa9euBWDMmDEMGjTIrxDDojb1lsQWLX3Dz7aqfiIi\nkRRzSVRBQQFHjhwhLS2t3Pa0tDTy8vIA6NixIyeeeKIf4YVNbeotiS1a+oafbVX9REQiKVauzquT\nu+66y+8QRKKS+oaISMOJuZGoNm3akJycTH5+frnt+fn5tG/f3qeowi9R6y21Fy1txM84ouUzEJHE\nEHNJVEpKCn379mXVqlWl25xzrFq1iv79+/sYWXglar2l9qKljfgZR7R8BiKSGGLydN7EiRMZM2YM\nffv2pV+/fsyaNYuioiLGjBkDwOTJk/nkk0949NFH/Q20gdVU73hWWFjIe++9h3MOgJ07d5Kbm0vr\n1q054YQTfI4uekRL3/Czrfp5bL/bqd/HF0k0MZlEjRw5koKCAqZPn05+fj69e/dmxYoVtG3bFoC8\nvDz27Nnjc5QNr6Z6x7P169czaNAgzAwzY9KkSQCMHj2ahQsX+hxd9IiWvuFnW/Xz2H63U7+PL5Jo\nYjKJAsjMzCQzM7PS1xYtWhThaCKnunrHs/POOy8h7/8XimjpG362Vb+O7Xc79fv4Iokm5uZEiYiI\niESDmBqJysrKYuvWrXXe74033qjX/rEsUeu+a9cuv0OImFjtF34eP5Hrnkh9QyTcrGQCYjQzs3OS\nkpLW1meYOikpKWGHuRO17sF6Zzjn1vkdSzjEQ7/w8/iqe/z2DZFIiZWRqMOBQIDFixfTs2fPOu+c\nlZXFtGnTQt4/liVq3bdu3cqVV14JcNjvWMIopvuFn8dP5LonSN8QiYhYSaIA6NmzJ3369KnzfiXD\n5aHuH8sSue6JIlb7hZ/HT+S6i0jDicmJ5dnZ2VxyySV06tSJpKQknn/+eb9Dioi7776bfv360bJl\nS9LS0vjxj3/Mjh07/A4rYubMmcNJJ51Eamoq55xzDu+8847fIUWVaOsXfn1f0fA5+NlW1U9EIicm\nk6jCwkJ69+7N3LlzMTO/w4mY7Oxsrr/+et566y1efvlliouLueCCCzhw4IDfoYXdU089xaRJk5gx\nYwYbN24kPT2dYcOGUVBQ4HdoUSOa+oWf35ffn4OfdVc/EYkw51zUP4A+gMvJyXEVmZlbtmzZUdvL\nWrx4satq/1i2b98+Z2YuOzu7yjLxUvezzz7b3XDDDaXPA4GA69Spk5s5c2al5XNychzggD4uCtpw\nOMFNdwcAABQeSURBVB7R3C9q831Fom1W9TnEa91rc+xE6Bt66BGpR0yORNXktttu46STTvI7jLD7\n8ssvMTNat27tdyhhVVxcTE5ODkOGDCndZmYMHTqUtWvX+hhZbIlUv0jk78vPuify5y7il5iaWF5b\nbdu2pVu3bn6HEVbOOW688UbOPfdcevXq5Xc4YVVQUMCRI0dIS0srtz0tLY3t27f7FFXsiVS/SOTv\ny8+6J/LnXhkz6wy0qcdbNCY2rmBUnA2rwDn3YW0Lx2USNW7cOMaNG+d3GGGVmZnJli1bShftE6lJ\nIvQLEfASqKSkpO2BQKBpqO/h9zpitaU4G1ZSUtJBM+te20QqLpOoeDd+/HiysrLIzs6mQ4cOfocT\ndm3atCE5OZn8/Pxy2/Pz82nfvr1PUUlVEvn78rPuify5V6JNIBBoGqtrqNWW4mxYwTXUmuKNYCqJ\nikfjx49n2bJlvPbaa3Tu3NnvcCIiJSWFvn37smrVKi655BLAO525atUqbrjhBp+jk4oS+fvys+6J\n/LlXJVbXUKstxem/mJxYXlhYSG5uLps2bQJg586d5ObmsmfPHsBbJ2Xo0KF+hhgWmZmZPP744zzx\nxBM0b96c/Px88vPzOXjwoN+hhd3EiRNZsGABjz32GNu2bWPs2LEUFRUxZswYv0OLGtHUL/z8vmr6\nHMLNz7qrn1QtVtbPioZ1zuoi2j/X+fPnk56eTqtWrWjVqhX9+/fnxRdfbLD3j8mRqPXr1zNo0CDM\nDDNj0qRJAIwePZqFCxdSUFDAzp07fY6y4c2fPx8zY+DAgeW2L1q0iFGjRvkTVISMHDmSgoICpk+f\nTn5+Pr1792bFihW0bdvW79CiRjT1Cz+/r5o+h3Dzs+7qJ5UrWT/rwQcfpF+/fsyaNYthw4axY8cO\n2rSpz9zzhleyztlVV13FZZdd5nc41YqFz/WEE05g5syZnHLKKTjneOSRRxg+fDibNm1qmFOLfq+x\nUJsH1ayHUxvxslZSKBK17omwFk6s9ws/j5/IdU/EvlHXdeb8bh8lalrvze84a/u5+h1nRa1bt3YL\nFy48ansofSMmT+eJiIjURm3WzxozZgyDBg3yK8SYFIvrkgUCAZ588kmKiorIyMhokPeMqdN5WVlZ\npRPU6qJkGYBQ949liVr3Xbt2+R1CxMRqv/Dz+Ilc90TqG1C79bM6duzoR2gxLZbWJdu8eTMZGRkc\nPHiQFi1asHTpUnr06NEg7x0rSVTjpKQkpk2bFvIb1Hf/WJaodQ+uS9LY7zjCKOb7hZ/HT/S6x3nf\nqJO77rrL7xAkjHr06EFubi779+9nyZIljBo1ijVr1jRIIhUrSdThQCAQ8hoTsbJGRTgkat2D631A\nbKyQG6qY7hd+Hj+R654gfaOU1s8Kj1j6XBs1akTXrl0BOOOMM3j77be5//77mTdvXv3fu97vEEHx\nvuZHOCRy3RNFrPYLP4+fyHVPNFo/Kzxi+XMNBAIcOnSoQd4r5iaW33333fTr14+WLVuSlpbGj3/8\nY3bs2OF3WBER7vUuol20r0fit2jpG36vc+Pn8f38DqLl+49GNa2fNXnyZEaPHu1vkEF+r3NWF7Gw\nLtmUKVPIzs7mgw8+YPPmzUyePJnXXnutZDS23mIuicrOzub666/nrbfe4uWXX6a4uJgLLriAAwcO\n+B1a2JWsd7FhwwZycnIYPHgww4cPT4gJ4yXrkcyYMYONGzeSnp7OsGHD/n979xMTxf3GcfwzW/jF\nhigJ/is2cDBNizYBkdgWQyJQ0mIaY0+eDOB1oYbWiyUptrHRkpo2plWIIaBJrbQxQT3QlGAT2ENb\n3bVaja2JsX+0RZSD1j/FGmd+B8umqAiOO/ud2Xm/kj2w7vI8X2ce9tnZ2Wc0OjpqOjXf8EttjM+5\n2blzpyzLSmts0/FNbgO/bH8/WrNmjbZt26bW1laVlpbqxx9/nDA/6+LFi75pUuLxuEpLS1VWVpac\nc7Z06VJt2rTJdGr3mer/1Q8uXbqk+vp6FRUVqaamRolEQv39/aqurk5NgOnOQjB500Pm4Vy+fNmx\nLMuJxWKTzoTw24yKVJps3sW4TFn7o855CeMsnHtNVRvp2DceNufGZPx01cVk2yAd8SeLTW1MLSh/\nN8kztUI5J+rKlSuyLEt5eXnJ+9atW5fxMz+8mHfhV0GcR+IHYa0NP3nQNghDbCAsAnVi+b0cx1Fz\nc7MqKiq0ePHi5P35+fnj70YyjpfzLvwqSPNI/CKMteE3k22DTI/tJ0GdoTZd5JlabmaoBbqJikaj\nOn36dHIDjcvkmR9ezrtA5ghjbfjNZNsg02P7ROBnqE0XeabWo85QC2wT1dTUpL6+PsViMeXn55tO\nJ228nHfhV0GaR+IHYa0NPzG5Ddj+kgI+Q226yDO13MxQC2QT1dTUpIMHD2pwcFCFhYWm0zEqlfMu\n/CrI80jSjdowz+Q2YPtPND6HKxaL6cMPP1QikdDw8LAOHDiQ/FvyICbneD1Krn6ZN/bBBx+opaVF\nzc3N+uijj+77d5N5bt26Vb29vfr555/15JNPavny5Wpra9Ozzz6bkt8fuBPLo9Go9u7dq88//1w5\nOTkaGRnRyMiIxsbGko9paWnxzcyPVPJ63oWfBWEeiWl+qQ3Tc25Mxp/ONsjE2H5neuzGowhSrpJ0\n9OhR7dq1SyUlJaZTeSCvR38E7khUR0eHLMtSZWXlhPu7u7tVV1cnSRoeHvbNzI9UGp93MTw8rNzc\nXBUXF6d23oWPrVmzRqOjo2ptbdXIyIiWLFniu3kkpvmlNuLxuKqqqmRZVnLOjSTV19erq6vL09im\n409nG2RibL+rra1VbW2tJPn+ixVByvX69etau3atOjs7tXnzZtPpPFBfX9+En3fv3q158+YpkUio\noqLisX9/4Joo27anfEx3d3caMkm/zs5O0ykYFY1GFY1GTafhW36pjRUrVkwrl0yMb3LdJmMH3bvv\nvqs9e/a4+nZWmDU2NmrVqlWqrq72bRN1r1SP/ghcEwUAQCrNnTtXzzzzjOk0AqWnp0fHjx9XPB43\nncq0eTH6I1BNVKbP/PBCWNcepneUQa0Lk/HDvPYw1cZ0NTY2qrGx0XQagXHhwgU1NzdrYGBA2dnZ\nptOZNk9Gf0x3tLnJm6SXIpHI+Dh2V7fHfX6Qb2Fd+7/rfsn0/ktd+DM+a8/o2pj0si8PuxTROL9c\npmSqXE3leeDAAScSiTjZ2dlOVlaWk5WV5ViWlbzPtm1f5PlfjY2NTmFhofPbb79N+hg3l30JypGo\nUMz88EJY1+5m3kcABbouTMYP89pDUhvwUE1NjU6ePDnhvoaGBi1atEgbN2703bcKvRz9EZQmStLd\nGRNHjhxRe3u7fv31V0nS888/r9bW1uS3GR7EL7M0HtejzjqRMmPtXs/5CLrxbdvR0fFIteH1vrFj\nxw5t27ZNFy9eVElJiT755BMtW7YsLfFNxZ7uvupFfOpkcjdu3NDZs2fHj1Ilx17k5eWpoKBAO3bs\nUG9vrwYGBgxnOnWufpCTk3PfOUU5OTmaPXu2796sR6NR7du3T4cOHUqO/pCk3NxczZgx47F/f+Dm\nRBUUFKitrU3Hjh1TIpFQdXW1Vq9eHYrzfYI2PyRVvJ7zkSn8VBtffPGFNmzYoPfee08//PCDSkpK\n9Oqrr2p0dDSjY5vcV6mTycXjcZWWlqqsrCw59mLp0qXatGmTpLvX5zx37pzhLO+aKle/8utrUkdH\nh/766y9VVlZqwYIFyduXX36ZmgDT/dzP5E0P+XzbcRwnLy/P6erqmvRzTj98Hptq0/lc33Eyc+2X\nL192LMtyYrHYpI9x89l20G5T1YXjPLw2vNw3XnzxRWf9+vXJn23bdp5++mmnra3N8/gmY99rsn01\nHfEni01tTC0ofzfJM7Xc1EbgjkT9l23b6unp0c2bN7V8+XJJdz+XraqqMpwZvJTqOR+ZyGRt3L59\nW4lEQi+//HLyPsuyVFNTo2+//TZjYz+IyX2VOgG8F6hzosadOnVK5eXlGhsb08yZM9Xb26vnnntO\nkrRgwQLD2cFLjpP6OR+ZxA+1MTo6qjt37mj+/PkT7p8/f77OnDmTsbHvZXJfpU6A9AhkE1VUVKQT\nJ07o6tWr2r9/v+rq6jQ0NKSioiJt2bLFdHrwkCdzPjIIteEfJvdV6uQut+cDjs/S8vu5tuSZWm7y\nC2QTlZWVpYULF0qSSktLdeTIEW3fvl3t7e2GM4OXmpqa1NfXp1gspvz8fNPp+JIfamPOnDl64okn\nkt+CGTcyMqKnnnoqY2P/l8l9lTqRJI1GIpGxtWvXuv76VSQSCcTF3ckztSKRyJht29P+Fkogm6h7\n2batW7dumU4DHvJyzkcmM1Eb2dnZKisr0+HDh5MjOBzH0eHDh7V+/fqMjT3O5L5KndzlOM7vlmU9\nJ2mO299h2/b/FIBZWuSZWrZtjzqO8/t0Hx+4JqqlpUUrV65UYWGhrl27pr1792pwcFD9/f2SpLff\nflt//vmn9uzZYzjT1AvC/BAveD3nI1P4qTbeeustNTQ0qKysTC+88II+/vhj3bx5Uw0NDRkd2+S+\nSp1M9O8L4bRfDAE3AtdEXbp0SfX19RoeHlZubq6Ki4vV39+v6upqSdLFixd1/vx5w1l6Ix6Pq6qq\nSpZlJeeHSFJ9fb26uroMZ+edjo4OWZalysrKCfd3d3errq7OTFI+5KfaWLNmjUZHR9Xa2qqRkREt\nWbJEX3/9tebOnZvRsU3uq9QJkH6Ba6I6Ozsf+u/d3d1pyiT9VqxYIdu2TaeRdmFcsxt+q41oNKpo\nNJrWmKZjm9xXqRMg/QI9JwoAAMAUmigAAAAXAvVxXl9fn6s5DuOzUtw+P8jCuvbxuSRhENS6MBk/\nzGsPU20AXrPGv+nlZ5ZlvRSJRL59nM/8I5FIaM8ZCOva/113ueM435nOxQuZUBcm47P2zK0NIF2C\nciTqH9u29dlnn2nRokWP/OS+vj698847rp8fZGFd+08//TQ+2M33c0keQ6DrwmT8MK89JLUBpMd0\nr1Rs8qaHXJF769atjmVZzptvvjnplZmDcgXpqWzZssVZtmyZM3PmTGfevHnO66+/7pw5c+ahz8mU\ntbe3tzvFxcXOrFmznFmzZjnl5eXOV199Nenjw36lepN1MTQ05KxatcpZsGCBY1mWc/DgQd/FD/Pa\nw1Ab3Lil6xboE8uPHj2qXbt2qaSkxHQqaRGLxfTGG2/o+++/18DAgG7fvq1XXnlFf//9t+nUPFdQ\nUKC2tjYdO3ZMiURC1dXVWr16dajO85ou03Vx48YNLVmyRDt37pRlWaGKH+a1A2EUlI/z7nP9+nWt\nXbtWnZ2d2rx5s+l00qKvr2/Cz7t379a8efOUSCRUUVFhKKv0eO211yb8/P7776u9vV3fffddqD6m\nnIof6qK2tla1tbWS7h7pDlP8MK8dCKPAHolqbGzUqlWrktOY/2vdunWqqqoykFV6XblyRZZlKS8v\nz3QqaWXbtnp6enTz5k2Vl5ebTsdXqAsASJ9AHonq6enR8ePHFY/HH/jv+fn5Gf8uzHEcNTc3q6Ki\nQosXLzadTlqcOnVK5eXlGhsb08yZM9Xb26uioiLTafkGdQEA6RW4JurChQtqbm7WwMCAsrOzH/iY\nLVu2pDmr9ItGozp9+nRy3kwYFBUV6cSJE7p69ar279+vuro6DQ0N0UiJugAAEwLXRCUSCV2+fFlL\nly5Nvqu+c+eOhoaG9Omnn+rWrVsZf0JlU1OT+vr6FIvFlJ+fbzqdtMnKytLChQslSaWlpTpy5Ii2\nb9+u9vZ2w5mZR10AQPoFromqqanRyZMnJ9zX0NCgRYsWaePGjRn/QtHU1KSDBw9qcHBQhYWFptMx\nyrZt3bp1y3QavhD2ugAAEwLXROXk5Nx3DlBOTo5mz56d/JZWS0uL/vjjD+3Zs8dEip6JRqPat2+f\nDh06pJycHI2MjEiScnNzNWPGDMPZeaulpUUrV65UYWGhrl27pr1792pwcFD9/f2mU/MFP9XFjRs3\ndPbs2eQRsXPnzunEiRPKy8tTQUGBp7FNxw/z2oEwClwT9SD3vsseHh7W+fPnDWXjnY6ODlmWpcrK\nygn3d3d3q66uzkxSaXLp0iXV19dreHhYubm5Ki4uVn9//wO/hYa7TNVFPB5XVVWVLMuSZVnasGGD\nJKm+vl5dXV0ZHT/MawfCKCOaqG+++WbCz93d3YYy8VYYr383rrOz03QKgWOqLlasWGF0XzUZP8xr\nB8IosHOiAAAATArUkSi3l/j45ZdfHuv5QRbWtYdpvUGtC5Pxw7z2MNUG4DUrCMP3LMsqjEQiZ2zb\ndn32dCQSCe1h7rCuPRKJjNm2/ZzjOL+bzsULmVAXJuOHfO0ZXRtAugSiiZLuvmBImvMYv+J/kv5J\nUTpBE9a1j2b6i0QG1IXJ+GFee8bXBpAOgWmiAAAA/IQTywEAAFygiQIAAHCBJgoAAMAFmigAAAAX\naKIAAABcoIkCAABwgSYKAADABZooAAAAF2iiAAAAXKCJAgAAcIEmCgAAwAWaKAAAABdoogAAAFyg\niQIAAHCBJgoAAMAFmigAAAAXaKIAAABcoIkCAABwgSYKAADABZooAAAAF2iiAAAAXKCJAgAAcIEm\nCgAAwAWaKAAAABdoogAAAFygiQIAAHCBJgoAAMAFmigAAAAXaKIAAABcoIkCAABwgSYKAADABZoo\nAAAAF2iiAAAAXKCJAgAAcIEmCgAAwAWaKAAAABdoogAAAFygiQIAAHCBJgoAAMAFmigAAAAXaKIA\nAABcoIkCAABwgSYKAADABZooAAAAF2iiAAAAXKCJAgAAcIEmCgAAwAWaKAAAABdoogAAAFygiQIA\nAHCBJgoAAMAFmigAAAAXaKIAAABcoIkCAABw4f/z3KJz6NWqHAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f4c676563c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "axes[0, 0].set_title"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3、任意 $p$ Lagrange 有限元空间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 三角形单元上的重心坐标"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAGHCAYAAABS9T4EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAEbBJREFUeJzt3WuM5Xddx/HPd3eb7U2lUCmlSotCociDCooLMc0GjFKC\nRCUYFQrFCmqJEiJpkYCADzAaUNIYBRtkJRWUS1CpNhoiFUgpEbACpWgruCCtlCIgYLtsd38+OLPt\n7O7MzpmZc/lfXq9ksrv/6283M+d9/rez1VoLAOO2Y9kDAGD5xAAAMQBADACIGAAQMQAgYgBAxACA\niAEAEQMAIgYARAwYgKq6tKoOV9XDlj2WpHvjgWmIAZ1SVU+sqldV1XduYrW28tUVXRsPbEgM6Jon\nJfmtJA/YxDpvTXJKa+3z8xkSDJ8Y0DU19YJVpyZJm/j2/IYEwycGdEZVvSrJ76388T9XzrsfqqqH\nVdWrV/58QVW9rar+J8kHV9Y76hz9yvJ/VFWfqar/q6q7quodVXXuGvs8st3vr6p9VfXVqvpaVf1p\nVZ18zLJ7q+qjVXV3Vd1aVS88sv4Uf7eHrmzzv6vqnqr6VFU9f9v/aDAju5Y9AFjl3UnOT/JzSV6c\n5Csr07+c+8/BvzPJvyf5zdx/FHHsOfofTrInyduT/FeS85JcnuT9VfWY1to9q5Y9st47knw2ycuS\nPC7JLyX50sp+UlU/mOS6JLcneWUmPzuvTHJXNrg+UFUPTvKRJIeSXLWyzsVJ3lxV39Fau+qE/yqw\nAGJAZ7TWPlVVH88kBn+9+hpA1X1nj/6ltXbJBpu6trX27tUTquq9SW5M8swkf77GOh9rrb1w1fJn\nJrksKzFI8pok9yZ5UmvtSyvLvCPJZ6b4q702k3Bd2Fr72sq0P6mqtyV5dVW9qbV2YL2Vq+pJmVwX\n+b5jZrUkl7fW3jTFGOCEnCaiT1qSDV/4Vr+wVtWuqnpgJu/6v5bJu/5ptvvBJA+qqtOrakeSpyT5\nqyMhWNnPZzM5WtjIzyR5b5KdVfWgI19J/iHJd60zpiPjf3CS30jyrEyOcK5Jcu7K788TAmal00cG\nVfVDSS5NsjeTb/6vZPLu7hWttVuXNjCW6XMbLbByrv/lmXzvnJOjTyd91zqrHXsn0ldXfj0jyelJ\nTkly2xrrrTVt9Vi+O5M7o16Y5JfXWKQlefAJNvH4JL/YWvv6qgvmXzjRPmErOh2DJFdmcqvhO5N8\nIslDkvxako9X1Y+01j69zMGxFHdPscwfJnlekj/I5M3D1zN50f3LrH80fGid6VPf3bSOI/u7Jsmf\nrbPMJ9ZbubW2+sjjyUn+bZvjgTV1PQavT/LzrbV7j0xYOU/7yUwu9D13WQNjbmbxsNYzk+xrrV1x\nZEJV7c7mnl1Y7c4k9yR5xBrzHrnBul9O8o0kO1tr/7jF/R/x1CR/v81twJo6fc2gtXbj6hCsTLst\nyc1JLljOqJizb638utUX7mTyLv/Y7+1fT7JzKxtrrR1O8r4kP1VVDzkyvaoekckL9EbrvjvJM6vq\nB46dv3KhelpPTfKvm1geptb1I4P1nJXkU8seBHPxsUxOzby2qv4iycFMLr5uxrVJLqmq/03y6SRP\nzOQC8F3bGNerk/x4khuq6o8z+dl5USZHqRdusO7LMrnu9ZGqunplTA/M5HrAk5NsGISqOj/Jgzxl\nzbz0LgZV9ZxMLgq+YtljYfZaax+tqlck+ZUkP5HJO/yHb3IzL87kNtBfSHJykg8l+bFMTrFs6TRU\na+3jVfXUJK9L8ttJvpDJcwaPSfLoDda9s6qekMnHbPx0kl/N5GaIm5NccaJ1V3l0kr/bythhGtVa\nfz5Pq6oenckFwU8muaj1afAMUlW9J8ljWmuPWvZYYDs6fc1gtao6K8nfZnLL37OEgEVb4+MpHpnk\naUnev5wRwez04shg5eOM/ynJ9yT50daa2+tYuKq6Pcm+TB5gOy+TU1knJXlca+0/ljcy2L7OXzNY\nuSXw2kxu63uKELBE12XyURkPSXIgyQ1JXi4EDEGnjwxWPgbgPZncUveM1pp7rAHmoOtHBr+f5CeT\n/E2SM6vq2atnttbW+sAxADaps0cG9Zpq2Zdkf47cDHjcfyXYWtvSQ0QAHK3bdxNdmsnHje3OoSS/\n21rbufpruYMDGI5ux+BgJo/ZHEglefOSRwMwWN09TXRetR137sjhuw+nqm5prT2+tTbNJ1YCsEmd\nPTI481tntitefEWuueaa7N69+1E7d+68tqpOWfa4AIaos0cGF1100b0f+MAHdibJ9ddfn4svvvjw\nwYMHrz906NDTHSEAzFZnjwxW27t3b6677rodJ5100l5HCACz14sYJIIAME+9iUEiCADz0qsYJIIA\nMA+9i0EiCACz1ssYJIIAMEu9jUEiCACz0usYJIIAMAu9j0EiCADbNYgYJIIAsB2DiUEiCABbNagY\nJIIAsBWDi0EiCACbNcgYJIIAsBmDjUEiCADTGnQMEkEAmMbgY5AIAsBGRhGDRBAATmQ0MUgEAWA9\no4pBIggAaxldDBJBADjWKGOQCALAaqONQSIIAEeMOgaJIAAkYpBEEADEYIUgAGMmBqsIAjBWYnAM\nQQDGSAzWIAjA2IjBOgQBGBMxOAFBAMZCDDYgCMAYiMEUBAEYOjGYkiAAQyYGmyAIwFCJwSYJAjBE\nYrAFggAMjRhskSAAQyIG2yAIwFCIwTYJAjAEYjADggD0nRjMiCAAfSYGMyQIQF+JwYwJAtBHYjAH\nggD0jRjMiSAAfSIGcyQIQF+IwZwJAtAHYrAAggB0nRgsiCAAXSYGCyQIQFeJwYIJAtBFYrAEggB0\njRgsiSAAXSIGSyQIQFeIwZIJAtAFYtABggAsmxh0hCAAyyQGHSIIwLKIQccIArAMYtBBggAsmhh0\nlCAAiyQGHSYIwKKIQccJArAIYtADggDMmxj0hCAA8yQGPSIIwLyIQc8IAjAPYtBDggDMmhj0lCAA\nsyQGPSYIwKyIQc8JAjALYjAAggBslxgMhCAA2yEGAyIIwFaJwcAIArAVYjBAggBslhgMlCAAmyEG\nAyYIwLTEYOAEAZiGGIyAIAAbEYOREATgRMRgRAQBWI8YjIwgAGsRgxESBOBYYjBSggCsJgYjJgjA\nEWIwcoIAJGJABAEQA1YIAoybGHAfQYDxEgOOIggwTmLAcQQBxkcMWJMgwLiIAesSBBgPMeCEBAHG\nQQzYkCDA8IkBUxEEGDYxYGqCAMMlBmyKIMAwiQGbJggwPGLAlggCDIsYsGWCAMMhBmyLIMAwiAHb\nJgjQf2LATAgC9JsYMDOCAP0lBsyUIEA/iQEzJwjQP2LAXAgC9IsYMDeCAP0hBsyVIEA/iAFzJwjQ\nfWLAQggCdJsYsDCCAN0lBiyUIEA3iQELJwjQPWLAUggCdIsYsDSCAN0hBiyVIEA3iAFLJwiwfGJA\nJwgCLJcY0BmCAMsjBnSKIMByiAGdIwiweGJAJwkCLJYY0FmCAIsjBnSaIMBiiAGdJwgwf2JALwgC\nzJcY0BuCAPMjBvSKIMB8iAG9Iwgwe2JALwkCzJYY0FuCALMjBvSaIMBsiAG9JwiwfWLAIAgCbI8Y\nMBiCAFsnBgyKIMDWiAGDIwiweWLAIAkCbI4YMFiCANMTAwZNEGA6YsDgCQJsTAwYBUGAExMDRkMQ\nYH1iwKgIAqxNDBgdQYDjiQGjJAhwNDFgtAQB7icGjJogwIQYMHqCAGIASQQBxABWCAJjJgawiiAw\nVmIAxxAExkgMYA2CwNiIAaxDEBgTMYATEATGQgxgA4LAGIgBTEEQGDoxgCkJAkMmBrAJgsBQiQFs\nkiAwRGIAWyAIDI0YwBYJAkMiBrANgsBQiAFskyAwBGIAMyAI9J0YwIwIAn0mBjBDgkBfiQHMmCDQ\nR2IAcyAI9I0YwJwIAn0iBjBHgkBfiAHMmSDQB2IACyAIdJ0YwIIIAl0mBrBAgkBXiQEsmCDQRWIA\nSyAIdI0YwJIIAl0iBrBEgkBXiAEsmSDQBWIAHSAILJsYQEcIAsskBtAhgsCyiAF0jCCwDGIAHSQI\nLJoYQEcJAoskBtBhgsCiiAF0nCCwCGIAPSAIzJsYQE8IAvMkBtAjgsC8iAH0jCAwD2IAPSQIzJoY\nQE8JArMkBtBjgsCsiAH0nCAwC2IAAyAIbJcYwEAIAtshBjAggsBWiQEMjCCwFWIAAyQIbJYYwEAJ\nApshBjBggsC0xAAGThCYhhjACAgCGxEDGAlB4ETEAEZEEFiPGMDICAJrEQMYIUHgWGIAIyUIrCYG\nMGKCwBFiACMnCCRiAEQQEANghSCMmxgA9xGE8RID4CiCME5iABxHEMZHDIA1CcK4iAGwLkEYDzEA\nTkgQxkEMgA0JwvCJATAVQRg2MQCmJgjDJQbApgjCMIkBsGmCMDxiAGyJIAyLGABbJgjDIQbAtgjC\nMIgBsG2C0H9iAMyEIPSbGAAzIwj9JQbATAlCP4kBMHOC0D9iAMyFIPSLGABzIwj9IQbAXAlCP4gB\nMHeC0H1iACyEIHSbGAALIwjdJQbAQglCN4kBsHCC0D1iACyFIHSLGABLIwjdIQbAUglCN4gBsHSC\nsHxiAHSCICyXGACdIQjLIwZApwjCcogB0DmCsHhiAHSSICyWGACdJQiLIwZApwnCYogB0HmCMH9i\nAPSCIMyXGAC9IQjzIwZArwjCfIgB0DuCMHtiAPSSIMyWGAC9JQizIwZArwnCbIgB0HuCsH1iAAyC\nIGyPGACDIQhbJwbAoAjC1ogBMDiCsHliAAySIGyOGACDJQjTEwNg0ARhOmIADJ4gbEwMgFEQhBMT\nA2A0BGF9YgCMiiCsTQyA0RGE44kBMEqCcDQxAEZLEO4nBsCoCcKEGACjJwhiAJBEEMQAYMWYgyAG\nAKuMNQhiAHCMMQZBDADWMLYgiAHAOsYUBDEAOIGxBEEMADYwhiCIAcAUhh4EMQCY0pCDIAYAmzDU\nIIgBwCYNMQhiALAFQwuCGABs0ZCCIAYA2zCUIIgBwDYNIQhiADADfQ+CGADMSJ+DIAYAM9TXIIgB\nwIz1MQhiADAHfQuCGADMSZ+CIAYAc9SXIIgBwJz1IQhiALAAXQ+CGAAsSJeDIAYAC9TVIIgBwIJ1\nMQhiALAEXQtC52Pw7W9/O1deeWXOOeecnHrqqdmzZ0/e9773LXtYANs2iyBU1WlV9Zqquq6qvlJV\nh6vquZvdTudj8LznPS9veMMbcskll+Sqq67Krl278rSnPS033HDDsocGsG0zCMKZSV6Z5NFJbkrS\ntjKOam1L683dRRdddO/rXve6nXv27MnrX//6vOQlL0mSHDhwII997GNz1lln5UMf+tCSRwkwG9df\nf30uvvjiwwcPHrz+0KFDT2+t3T3NelV1UpIzWmt3VtXjk/xzkktba2/dzP47fWTwrne9K7t27coL\nXvCC+6bt3r07l112WT784Q/ni1/84hJHBzA7Wz1CaK0dbK3dud39dzoGN910U84///ycfvrpR01/\nwhOecN98gKFY5kXlTsfgjjvuyNlnn33c9LPPPjuttdx+++1LGBXA/CwrCJ2Owd13353du3cfN/3k\nk0++bz7A0CwjCJ2OwSmnnJIDBw4cN/2ee+65bz7AEC06CJ2NwS233LLjm9/8Zvbv33/cvDvuuCNJ\n8tCHPnTRwwJYmLWCUFWPrKrfqaq3z3JfnY3BXaffVfu/uD+33npr3vjGNx4178Ybb0xV5cILL1zS\n6AAWY3UQqupjST6T3Xlpzs3PznI/nY1BLk3asyfPQFx++eW57bbbkkyeSN63b1/27NmTc845Z4kD\nBFiMvXv35uqrr97RWrsgj8uOvDS78vzZvn7vmuXGZu68JBck7ZaW5zznObnsssuyb9++7N+/P295\ny1uWPTqAhbn55puz49QdOXzx4eSko+dV1YuSPCDJkXfIz6iq7135/VWttW9stP3OPoFcr6nJwO5N\nclWy45uTCJ522mk599xzD59xxhndHDjAHNxyyy077jr9rsql909rr2qVJFX1uSQPW2fVh7fWPr/R\n9jsbAwAWp7vXDABYGDEAQAwAEAMAIgYARAwAiBgAEDEAIGIAQMQAgIgBABEDACIGAEQMAIgYABAx\nACBiAEDEAICIAQARAwAiBgBEDACIGAAQMQAgYgBAxACAiAEAEQMAkvw/nbt+wtdQYjAAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd6381673c8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from fealpy.mesh.TriangleMesh import TriangleMesh\n",
    "\n",
    "#from matplotlib.font_manager import FontProperties\n",
    "#cfont = FontProperties('SimHei')\n",
    "\n",
    "point = np.array([\n",
    "    [0,0],\n",
    "    [1,0],\n",
    "    [0,1]], dtype=np.float)\n",
    "cell = np.array([[0, 1, 2]], dtype=np.int)\n",
    "\n",
    "mesh = TriangleMesh(point, cell)\n",
    "fig = plt.figure()\n",
    "axes = fig.gca()\n",
    "mesh.add_plot(axes, cellcolor='w')\n",
    "mesh.find_point(axes, showindex=True, fontsize=12, color='g', markersize=25)\n",
    "axes.set_title('triangle $\\\\tau$')\n",
    "plt.tight_layout(pad=1, w_pad=1, h_pad=1.0)\n",
    "plt.savefig('/home/why/ppt/figures/triangle.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定三角形单元 $\\tau$, 其三个顶点 $\\mathbf x_0 :=(x_0,y_0)$, $\\mathbf x_1 :=(x_1,y_1)$ 和 $\\mathbf x_2 :=(x_2,y_2)$ 逆时针排列, 且不在同一条直线上, 那么向量 $\\vec{\\mathbf x_0\\mathbf x_1}$ 和 $\\vec{\\mathbf x_0\\mathbf x_2}$ 是线性无关的. 这等价于矩阵\n",
    "\n",
    "$$\n",
    "A = \n",
    "\\begin{pmatrix}\n",
    "x_0 & x_1 & x_2 \\\\\n",
    "y_0 & y_1 & y_2 \\\\\n",
    "1   & 1   & 1 \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "非奇异. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "任给一点 $\\mathbf{x}:=(x,y)\\in\\tau$, 求解下面的线性方程组\n",
    "\n",
    "$$\n",
    "A \n",
    "\\begin{pmatrix}\n",
    "\\lambda_0 \\\\\n",
    "\\lambda_1\\\\\n",
    "\\lambda_2  \n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "x \\\\\n",
    "y\\\\\n",
    "1  \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "可得唯一的一组解$\\lambda_0,\\lambda_1,\\lambda_2$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因此对任意二维点 $\\mathbf{x}\\in\\tau$, 有\n",
    "\n",
    "$$\n",
    "\\mathbf{x}=\\lambda_0(\\mathbf{x})\\mathbf{x}_0 + \\lambda_1(\\mathbf{x})\\mathbf{x}_1 + \\lambda_2(\\mathbf{x})\\mathbf{x}_2 \n",
    "\\text{ 且 } \\lambda_0(\\mathbf{x}) + \\lambda_1(\\mathbf{x}) + \\lambda_2(\\mathbf{x}) = 1. \n",
    "$$\n",
    "\n",
    "$\\lambda_0,\\lambda_1,\\lambda_2$ 称为点 $\\mathbf{x}$ 关于点 $\\mathbf{x}_0,\\mathbf{x}_1$ 和$\\mathbf{x}_2$ 的**重心坐标**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "易知, $\\lambda_0, \\lambda_1, \\lambda_2$ 都是关于 $\\mathbf x$ 的线性函数, 且有\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\lambda_0(\\mathbf x_0) = 1,\\quad & \\lambda_0(\\mathbf x_1) = 0,\\quad& \\lambda_0(\\mathbf x_2) = 0\\\\\n",
    "\\lambda_1(\\mathbf x_0) = 0,\\quad & \\lambda_1(\\mathbf x_1) = 1,\\quad& \\lambda_1(\\mathbf x_2) = 0\\\\\n",
    "\\lambda_2(\\mathbf x_0) = 0,\\quad & \\lambda_2(\\mathbf x_1) = 0,\\quad & \\lambda_2(\\mathbf x_2) = 1\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\lambda_0, \\lambda_1, \\lambda_2$ 关于 $\\mathbf x$ 的梯度分别为:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla\\lambda_0 = \\frac{1}{2|\\tau|}(\\mathbf x_2 - \\mathbf x_1)W\\\\\n",
    "\\nabla\\lambda_1 = \\frac{1}{2|\\tau|}(\\mathbf x_0 - \\mathbf x_2)W\\\\\n",
    "\\nabla\\lambda_2 = \\frac{1}{2|\\tau|}(\\mathbf x_1 - \\mathbf x_0)W\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "其中 \n",
    "\n",
    "$$\n",
    "W = \\begin{pmatrix}\n",
    "0 & 1\\\\ -1 & 0 \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**目标:** \n",
    "\n",
    "给定有 $NC$ 个三角形单元的网格, 计算所有单元上的 $\\nabla \\lambda_0$, $\\nabla \\lambda_1$ 和 $\\nabla \\lambda_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里需要一个三维数组来存储所有单元上的梯度. \n",
    "\n",
    "记三维数组为 $Dlambda$, 则 $Dlambda[i, j, k]$ 存储第 $i$ 个单元上的第 $j$ 个重心坐标的梯度的第 $k$ 个分量的值, 其中 $0 \\leq i \\leq NC-1$, $ 0 \\leq j \\leq 2$, $ 0 \\leq k \\leq 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dlambda:\n",
      " [[[ 1. -1.]\n",
      "  [ 0.  1.]\n",
      "  [-1.  0.]]\n",
      "\n",
      " [[-1.  1.]\n",
      "  [ 0. -1.]\n",
      "  [ 1.  0.]]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "point = np.array(\n",
    "    [(0.0, 0.0),\n",
    "     (1.0, 0.0),\n",
    "     (1.0, 1.0), \n",
    "     (0.0, 1.0)], dtype=np.float)\n",
    "cell = np.array([\n",
    "        (1, 2, 0), \n",
    "        (3, 0, 2)], dtype=np.int)\n",
    "NC = cell.shape[0]\n",
    "v0 = point[cell[:, 2], :] - point[cell[:, 1], :] # x_2 - x_1\n",
    "v1 = point[cell[:, 0], :] - point[cell[:, 2], :] # x_0 - x_2\n",
    "v2 = point[cell[:, 1], :] - point[cell[:, 0], :] # x_1 - x_0\n",
    "nv = np.cross(v2, -v1)\n",
    "\n",
    "Dlambda = np.zeros((NC, 3, 2), dtype=np.float)\n",
    "length = nv\n",
    "W = np.array([[0, 1], [-1, 0]], dtype=np.int)\n",
    "Dlambda[:,0,:] = v0@W/length.reshape(-1, 1)\n",
    "Dlambda[:,1,:] = v1@W/length.reshape(-1, 1)\n",
    "Dlambda[:,2,:] = v2@W/length.reshape(-1, 1)\n",
    "print('Dlambda:\\n', Dlambda)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2  三角形单元上的 $p$ 次基函数公式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定三角形单元上的一个重心坐标 $(\\lambda_0, \\lambda_1, \\lambda_2)$, 所有 $p\\geq 1$ 次基函数的计算公式如下:\n",
    "\n",
    "$$\n",
    "\\phi_{m,n,k} = \\frac{p^p}{m!n!k!}\\prod_{l_0 = 0}^{m - 1}\n",
    "(\\lambda_0 - \\frac{l_0}{p}) \\prod_{l_1 = 0}^{n-1}(\\lambda_1 -\n",
    "\\frac{l_1}{p}) \\prod_{l_2=0}^{k-1}(\\lambda_2 - \\frac{l_2}{p}).\n",
    "$$\n",
    "\n",
    "其中 $ m\\geq 0$, $n\\geq 0$, $ k \\geq 0$, 且 $m+n+k=p$, 这里规定:\n",
    "\n",
    "$$\n",
    " \\prod_{l_i=0}^{-1}(\\lambda_i - \\frac{l_i}{p}) := 1,\\quad i=0, 1, 2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**注:**\n",
    "\n",
    "三维情形的基函数表示如下:\n",
    "$$\n",
    "\\phi_{m,n,k,v} = \\frac{p^p}{m!n!k!v!}\\prod_{l_0 = 0}^{m - 1}\n",
    "(\\lambda_0 - \\frac{l_0}{p}) \\prod_{l_1 = 0}^{n-1}(\\lambda_1 -\n",
    "\\frac{l_1}{p}) \\prod_{l_2=0}^{k-1}(\\lambda_2 - \\frac{l_2}{p})\\prod_{l_3=0}^{v-1}(\\lambda_3 - \\frac{l_3}{p}).\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 基函数的矩阵计算公式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**目标:**\n",
    "\n",
    "给定重心坐标 $(\\lambda_0, \\lambda_1, \\lambda_2)$, 利用**面向数组**的方式计算 $p$ 次元所有基函数在 $(\\lambda_0, \\lambda_1, \\lambda_2)$ 处的**函数值**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造向量：\n",
    "$$\n",
    "P = ( \\frac{1}{0!},  \\frac{1}{1!}, \\frac{1}{2!}, \\cdots, \\frac{1}{p!})\n",
    "$$\n",
    "\n",
    "构造矩阵：\n",
    "$$\n",
    "A :=                                                                            \n",
    "\\begin{pmatrix}  \n",
    "1  &  1  & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\                                             \n",
    "\\lambda_0 - \\frac{1}{p} & \\lambda_1 - \\frac{1}{p} & \\lambda_2 - \\frac{1}{p}\\\\   \n",
    "\\vdots & \\vdots & \\vdots \\\\                                                     \n",
    "\\lambda_0 - \\frac{p - 1}{p} & \\lambda_1 - \\frac{p - 1}{p} & \\lambda_2 - \\frac{p - 1}{p}\n",
    "\\end{pmatrix}                                                                   \n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对 $A$ 的每一列做累乘运算, 并左乘由 $P$ 形成的对角矩阵, 得矩阵:\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "易知, 只需从 $B$ 的每一列中各选择一项相乘(要求三项次数之和为 $p$), 再乘以 $p^p$ 即可得到相应的基函数, 其中取法共有 \n",
    "\n",
    "$$\n",
    "n_{dof} = {(p+1)(p+2)\\over 2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "即需要构造 $n_{dof}\\times 3$ 的指标矩阵 $I$, 则基函数可写成如下形式 \n",
    "\n",
    "$$\n",
    "\\phi_i = p^pB_{I_{i,0}, 0}B_{I_{i, 1},1}B_{I_{i, 2}, 2}, \\quad i = 0, 1, \\cdots, n_{dof}\n",
    "$$\n",
    "\n",
    "对应 `Python` 的代码如下:\n",
    "\n",
    "```Python\n",
    "phi = p**p*np.prod(B[I, [0, 1, 2]], axis=1)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlAAAAGICAYAAACOZ96aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVVX+//HX2kiBmBRekBgveGmgixA0hkYpaIr1Vcem\n/FaDilpZiKYy1WihOd6yfkU6qUwpaOGFmseYVmfUZEzJvJIy+Q0v5SVNPcbMaAqYyFm/PxBGFJUD\n55x9Dnyej8d6CPvss9fbczj6Ye2111Zaa4QQQgghRM0ZZgcQQgghhPA0UkAJIYQQQthJCighhBBC\nCDtJASWEEEIIYScpoIQQQggh7CQFlBBCCCGEnaSAEkIIIYSwkxRQQgghhBB2kgJKCCGEEMJOUkAJ\nIYQQQthJCighhBBCCDtJASWEEEIIYScpoIQQQggh7CQFlBBCCCGEnaSAEkIIIYSwkxRQQgghhBB2\nkgJKCCGEEMJOUkAJIYQQQthJCighhBBCCDtJASWEEEIIYScpoIQQQggh7CQFlBBCCCGEnaSAEkII\nIYSwkxRQQgghhBB2kgJKCCGEEMJOUkAJIYQQQthJCighhBBCCDtJASWEEEIIYScpoIQQQggh7CQF\nlBBCCCGEnaSAEkIIIYSwkxRQQgghhBB2amR2ACGEEMIMSqk2QHOzc9jhBuC82SHs4Gl5AQq11j/U\nZEcpoIQQQjQ4Sqk2hmHstdlsPmZnqSnDMLDZbGbHqDFPywtgGMY5pdSva1JESQElhBCiIWpus9l8\nsrKyCAsLMzvLdVksFlJTU5G8zlNQUEBCQoIP5aOSUkAJIYQQVxMWFkZkZKTZMa6roKAAkLzuRCaR\nCyGEEMLtxcbGMn78eLNjVJICSgghhBDCTlJACSGEEELYSQooIYQQopZiY2MZM2YM48aNIyAggFat\nWrFw4UKKi4sZPnw4TZs2pVOnTqxevdrsqFXExsby/PPP89JLL9GsWTOCgoKYMmWK2bHs8tlnn3Hz\nzTezbNkyU/qXAkoIIYSog/fff58WLVqwfft2xowZw7PPPstjjz3Gfffdx86dO+nduzdDhgzh3Llz\nZket4v3336dJkyZs27aN119/nT/96U/k5OSYHatGli5dyu9//3uWLVvGE088YUoGKaCEEEKIOggP\nD2fixIl06NCBP/7xj/j4+NCiRQtGjBhBhw4dmDRpEoWFhfzzn/80O2oVnTt3JjU1lQ4dOjB48GDu\nuecejyig5s2bR3JyMp9++il9+/Y1LYcsYyCEEELUQefOnSu/NgyDZs2acdddd1VuCwwMBODkyZMu\nz3Ytl+YGCAoKcruMl/voo4/46aef2LRpE1FRUaZmkREoIYQQog68vb2rfK+UumIb4HarcleX290y\nXi4yMpIWLVqwcOFCs6NIASWEEEIIz9ChQwfWr1/PypUrGT16tKlZ5BSeEEIIITxGx44dWb9+PbGx\nsTRq1Ii0tDRTckgBJYQQQtSSUqrW28zkbnlq4tLMt912Gzk5OZVF1BtvvOHyPFJACSGEELX0j3/8\n44ptBw4cuGJbWVmZK+LUWHW5V6xYYUKSmrs8c2hoKMePHzcpjcyBEkIIIYSwmxRQQgghhBB2kgJK\nCCGEEMJOMgdKCCGE21BKtQGau6CrUACLxUJBQYELuqubTZs2AZLXmQ4ePGjX/kpr7aQoQgghRM0p\npdoYhrHXZrP5uKI/wzDcfuHIS0le57uYuavWesv19pURKCGEEO6iuc1m88nKyiIsLMypHVksFlJT\nU3FFX44geZ2voKCAhIQEgPM12V8KKCGEEG4lLCyMyMhIp/ZRcVrJFX05guR1PzKJXAghhBDCTlJA\nCSGEcHtz584lJCQEX19foqOj2b59u9mRHO61117DMAzGjx9vdpRr8rT3wll5pYASQgjh1rKzs0lJ\nSWHKlCns3LmT8PBw+vTpQ2FhodnRHGb79u28++67hIeHmx3lmjztvXBmXimghBBCuLW0tDRGjhzJ\nkCFDCA0NJT09ncaNG5ORkWF2NIc4e/YsCQkJLFiwgJtvvtnsONfkae+FM/NKASWEEMJtlZaWkpeX\nR8+ePSu3KaXo1asXmzdvBiAxMZHY2FizItbZqFGj6NevH3FxcWZHuaaavBfuxNl55So8IYQQbquw\nsJCysjICAwOrbA8MDGTv3r0A3HrrrWZEc4jly5eza9cuduzYYXaU66rJe+FOnJ1XCighhBAebcaM\nGWZHqJWjR48yduxY1q1bh7e3t9lxhJ2kgBJCCOG2mjdvjpeXF1artcp2q9VKq1atTErlGHl5efz0\n009ERkZScVeQsrIyNm7cyDvvvMMvv/yCUsrklP/lae+Fs/PKHCghhBBuy9vbm6ioKHJyciq3aa3J\nycmhW7duJiaru169evHNN9+wa9cu8vPzyc/P55577iEhIYH8/Hy3Kp7A894LZ+eVESghhBBubfz4\n8SQmJhIVFUWXLl1IS0ujuLiYxMREACZMmMCxY8dYvHixuUHt5Ofnx+23337FtmbNmrnt7U+u9164\nG2fmlQJKCCGEWxs0aBCFhYVMmjQJq9VKREQEa9asoUWLFgCcOHGCI0eOmJzSMdxt1Oly13sv3I0z\n80oBJYQQwu0lJSWRlJRU7WOZmZkuTuM8//jHP8yOcF3Xei/ckbPyyhwoIYQQQgg7SQElhBBCCGEn\nKaCEEEIIIewkc6CEEEK4FYvFQkFBgVP72LRpk8v6cgTJ63wHDx60a39VsXiXEEIIYSalVLRhGJtt\nNptL+jMMA1f15QiS1/kuZu6qtd5yvX1lBEoIIYS7OG+z2cjKynL6OkgWi4XU1FSX9OUIktf5CgoK\nSEhIADhfk/2lgBJCCOFWwsLCiIyMdGofFaeVXNGXI0he9yOTyIUQQggh7CQFlBBCCLeWm5tL//79\nCQ4OxjAMVq1aZXYkh0hPTyc8PBx/f3/8/f3p1q0bq1evNjvWNXniezF37lxCQkLw9fUlOjqa7du3\nO+S4UkAJIYRwa0VFRURERDBv3jy3v9WJPVq3bs2sWbP4+uuvycvLIy4ujgEDBrj1VWue9l5kZ2eT\nkpLClClT2LlzJ+Hh4fTp04fCwsI6H1vmQAkhhHBr8fHxxMfHA1Cfrhx/+OGHq3w/bdo05s+fz5Yt\nW9x24rWnvRdpaWmMHDmSIUOGAOWjfp999hkZGRm8+OKLdTq2jEAJIYTwaK+++iohISFmx6gTm83G\n8uXLKS4upmvXrmbHqRdKS0vJy8ujZ8+elduUUvTq1YvNmzfX+fgyAiWEEMKjtWjRgo4dO5odo1Z2\n795N165dOXfuHDfddBMrVqwgNDTU7Fj1QmFhIWVlZQQGBlbZHhgYyN69e+t8fBmBEkII4dFGjRrF\n559/bnaMWgkNDSU/P59t27bx3HPPMWTIEPbs2WN2LFEDMgIlhBBCmKRRo0a0b98egLvvvptt27Yx\ne/Zs5s+fb3Iyz9e8eXO8vLywWq1VtlutVlq1alXn48sIlBBCCOEmbDYbv/zyi9kx6gVvb2+ioqLI\nycmp3Ka1Jicnh27dutX5+DICJYQQwq0VFRXx3XffVV71deDAAfLz8wkICKB169bMnTuXFStWsG7d\nOpOT2mfixIn07duXNm3acObMGZYsWcKGDRtYu3at2dGu6nrvhbsZP348iYmJREVF0aVLF9LS0igu\nLiYxMbHOx5YCSgghhFvbsWMHsbGxKKVQSpGSkgLA0KFDycjIoLCwkAMHDpic0n4nT55k6NChHD9+\nHH9/fzp37szatWuJi4szO9pVXe+9cDeDBg2isLCQSZMmYbVaiYiIYM2aNbRo0aLOx5YCSgghhFvr\n3r07Npvtqo9PnjyZyZMnuzCRYyxYsMDsCHa73nvhjpKSkkhKSnL4cWUOlBBCCCGEnaSAEkIIIYSw\nkxRQQgghhBB2kjlQQggh3IrFYnH6DXU3bdrksr4cQfI638GDB+3aX3nCzQCFEELUf0qpaMMwNrtq\nkrJhGB41IVryOt/FzF211luut6+MQAkhhHAX5202G1lZWYSFhTm1I4vFQmpqqkv6cgTJ63wFBQUk\nJCQAnK/J/lJACSGEcCthYWFERkY6tY+K00qu6MsRJK/7kUnkQgghhBB2kgJKCCGE25o5cyZdunSh\nadOmBAYGMnDgQPbt22d2LIdIT08nPDwcf39//P396datG6tXrzY71nXNnTuXkJAQfH19iY6OZvv2\n7WZHuqrc3Fz69+9PcHAwhmGwatUqhx1bCighhBBuKzc3l9GjR7N161bWrVtHaWkpvXv3pqSkxOxo\ndda6dWtmzZrF119/TV5eHnFxcQwYMMCtr1rLzs4mJSWFKVOmsHPnTsLDw+nTpw+FhYVmR6tWUVER\nERERzJs3D6WUQ48tc6CEEEK4LYvFUuX7RYsW0bJlS/Ly8oiJiTEplWM8/PDDVb6fNm0a8+fPZ8uW\nLW478TotLY2RI0cyZMgQoHwU7bPPPiMjI4MXX3zR5HRXio+PJz4+HgBHrzogI1BCCCE8xqlTp1BK\nERAQULlt2LBhxMbGmpiq7mw2G8uXL6e4uJiuXbuaHadapaWl5OXl0bNnz8ptSil69erF5s2bTUxm\nDhmBEkII4RG01owdO5aYmBhuv/32yu1BQUEOH11wld27d9O1a1fOnTvHTTfdxIoVKwgNDTU7VrUK\nCwspKysjMDCwyvbAwED27t1rUirzSAElhBDCIyQlJfHtt99WrnJdYcaMGSYlqrvQ0FDy8/M5ffo0\nf/3rXxkyZAgbN2502yJK/JcUUEIIIdxecnIyFouF3NxcgoKCzI7jMI0aNaJ9+/YA3H333Wzbto3Z\ns2czf/58k5NdqXnz5nh5eWG1Wqtst1qttGrVyqRU5pE5UEIIIdxacnIyK1euZP369bRp08bsOE5l\ns9n45ZdfzI5RLW9vb6KiosjJyancprUmJyeHbt26mZjMHDICJYQQwm0lJSWxbNkyVq1ahZ+fX+Xo\nh7+/Pz4+PgBMnDiRH3/8kcWLF5sZ1W4TJ06kb9++tGnThjNnzrBkyRI2bNjA2rVrzY52VePHjycx\nMZGoqCi6dOlCWloaxcXFJCYmmh2tWkVFRXz33XeVc+QOHDhAfn4+AQEBtG7duk7HlgJKCCGE20pP\nT0cpRY8ePapsz8zMrLyU/vjx4xw5csSEdHVz8uRJhg4dyvHjx/H396dz586sXbuWuLg4s6Nd1aBB\ngygsLGTSpElYrVYiIiJYs2YNLVq0MDtatXbs2EFsbCxKKZRSpKSkADB06FAyMjLqdGwpoIQQQrgt\nm8123X0yMzNdkMTxFixYYHaEWklKSiIpKcnsGDXSvXv3Gv0M1YbMgRJCCCGEsJNbF1BKqXuUUu8o\npXYrpc4qpQ4rpbKVUp3MziaEu1JK+Smlpiil/q6U+pdSyqaUGnKd5zRSSn17cd/xrsoqhBCeyq0L\nKOAlYCCwDhgD/AV4APhaKXX7tZ4oRAPWHEgFQoFdQE1WGBwDtK7hvkII0eC5+xyoN4EntNYXKjYo\npT4EvgH+CFzzt2ohGqhjQCut9UmlVBRwzVulK6VaUl5wvQZMdUE+Ia7JYrE4/Ya6FYtxuqIvR5C8\nznfw4EG79leeuPy9UmoHoLXWvzE7ixDu7JICKlFr/f5V9smgfLTqSeAA8Aet9VuuSylEOaVUtGEY\nm5016fdyhmE4bYKxM0he57uYuavWesv19nX3EairCQR2mx1CCE+nlOpC+UhuN+T0nTDfeZvNRlZW\nFmFhYU7tyGKxkJqa6pK+HEHyOl9BQQEJCQkA52uyv8cVUEqpBCAYeMXsLELUA38Glmmttyml2pod\nRgiAsLAwIiMjndpHxWklV/TlCJLX/XhUAaWUCgXeATYB1Z6OEELUjFJqGHAH5RdqCCGEsIO7X4VX\nSSkVCHwG/Ad4THvi5C0h3IRS6iZgBvC61vqY2XmEuJr09HTCw8Px9/fH39+fbt26sXr1arNjOcTM\nmTPp0qULTZs2JTAwkIEDB7Jv3z6zY12Vp+WtMHfuXEJCQvD19SU6Oprt2695XU2NeUQBpZRqCqwG\nmgLxWusTJkcSwtO9AHgDHyql2l48fVdxY6hbLm7zNi+eEOVat27NrFmz+Prrr8nLyyMuLo4BAwZ4\nzJVd15Kbm8vo0aPZunUr69ato7S0lN69e1NSUmJ2tGp5Wl6A7OxsUlJSmDJlCjt37iQ8PJw+ffpQ\nWFhY94Nrrd26ATcCG4EzQBez80iT5kkNiAJswJDLtmcCZRcfu7SVXfJnZ7PzS2tYDYgEdF5enr6W\ngIAAnZGRcc19ricrK0vXpC9X+umnn7RSSufm5l7xmOStnXvvvVePGTOm8nubzaaDg4P1rFmzrtg3\nLy9PU34xTaSuwc+rW49AKaUM4EPgXuBRrfU2kyMJUV/Mpnzu028vac8AivLi6reAfYuiCOFkNpuN\n5cuXU1xcTLdu3QBITEwkNjbW5GSOcerUKZRSBAQEmB2lRtw9b2lpKXl5efTs2bNym1KKXr16sXnz\n5jof390nkb8F9ANWAc2VUr+/9EGt9RJTUgnh5pRSo4CbKb9iFaC/UqriFN0crfUuylcpv/Q5FVfh\n/Z/W+hPXJBXi+nbv3k3Xrl05d+4cN910EytWrODXv/41ALfeeqvJ6RxDa83YsWOJiYnh9tvd/0Yb\nnpC3sLCQsrIyAgMDq2wPDAxk7969dT6+uxdQ4ZQPp/W72C4nBZQQ1fsD0Obi15ry0aaKq+0+oPyU\neHXk4gzhdkJDQ8nPz+f06dP89a9/ZciQIWzcuJHQ0FBmzJhhdjyHSEpK4ttvv61cwdvdeVpeZ3Db\nAkpNUZpXL36TSRmH+Uhr/YSZmYRwa0pVFj/doWwDfGjPZ0ZrfRjwcko2IeqgUaNGtG/fHoC7776b\nbdu2MXv2bObPn29yMsdITk7GYrGQm5tLUFCQ2XGuy1PyNm/eHC8vL6xWa5XtVquVVq1a1fn4bj0H\nCoBSoPyau0Om5hDCQ5QAO8u/PGRmDiGcxWaz8csvv5gdwyGSk5NZuXIl69evp02bNtd/gsk8Ka+3\ntzdRUVHk5ORUbtNak5OTUzmHri7cdgQKKC+eLMAvKGChyWmEcHslwGjgTPkvR/KZER5v4sSJ9O3b\nlzZt2nDmzBmWLFnChg0bWLt2LQATJkzg2LFjLF682OSk9ktKSmLZsmWsWrUKPz+/ypESf39/fHx8\nTE53JU/LCzB+/HgSExOJioqiS5cupKWlUVxcTGJiYp2P7b4FVCZ4/eRFWUkZSimb1vpXwHdmxxLC\nXXUHvvHy4rTNhtZaUX5z4D+ZHEuIOjl58iRDhw7l+PHj+Pv707lzZ9auXUtcXBwAJ06c4MiRIyan\nrJ309HSUUvTo0aPK9szMTIYMGWJOqGvwtLwAgwYNorCwkEmTJmG1WomIiGDNmjW0aNGizsd22wKq\neVFz/dSYp9Tvf/97xo4da6xfv361Uipea/2F2dmEcEcbgT+MG8fIkSPJzs7mlVdemaKUQmstRZTw\nWAsWLLjm45mZmS5K4ng2m83sCHbxtLwVkpKSSEpKcvhx3XYOVFhYmG3mzJnceeedfPLJJ0ZsbKy3\nYRirlVI9zM4mhLuaOnUqHTt25OWXX2batGkAU5RSk8zOJYQQ9Y3bFlCX8vX1lSJKCDtJESWEEM7j\nEQUUSBElRG1IESWEEM7htnOgqlNRRPXr189b5kQJUTMvv/wygMyJEh7DYrE4/WbBFQtAuqIvR5C8\nznfwoH13r1Jau+fCww888MCFjRs3VruoX0lJCf369bOtX7++1GazSRElBKCU0iUlJVe9nHj69Om8\n8sorAJOliBLuSCkVbRjGZldNVjYMw6MmRkte57uYuavWesv19vWoEagKMhIlhP1kJEp4gPM2m42s\nrCzCwsKc2pHFYiE1NdUlfTmC5HW+goICEhISAM7XZH+PLKBAiighakOKKOEJwsLCiIyMdGofFaeV\nXNGXI0he9+Mxk8irIxPLhbCfTCwXQoi68+gCCqSIEqI2pIgSnuq1117DMAzGjx9vdpQ6mzlzJl26\ndKFp06YEBgYycOBA9u3bZ3asa8rNzaV///4EBwdjGAarVq0yO9I1OTOvxxdQIEWUELUhRZTwNNu3\nb+fdd98lPDzc7CgOkZuby+jRo9m6dSvr1q2jtLSU3r17U1JSYna0qyoqKiIiIoJ58+ahlDI7znU5\nM6/HzoG6nMyJEsJ+MidKeIqzZ8+SkJDAggULmDp1qtlxHMJisVT5ftGiRbRs2ZK8vDxiYmJMSnVt\n8fHxxMfHA+CuV/Ffypl568UIVAUZiRLCfjISJTzBqFGj6NevX+VNhC81bNgwYmNjTUjlWKdOnUIp\nRUBAgNlRRA3UqwIKpIgSojakiBLubPny5ezatYuZM2dW+3hQUBBt27Z1cSrH0lozduxYYmJiuP32\n282OI2qg3pzCu5SczhPCfnI6T7ijo0ePMnbsWNatW4e3t3e1+8yYMcPFqRwvKSmJb7/9tnIFb+H+\n6t0IVAUZiRLCfjISJdxNXl4eP/30E5GRkXh7e+Pt7c2GDRuYPXs2N9xwg0fMw7me5ORkLBYLX3zx\nBUFBQWbHETVUbwsokCJKNEx5eXnEx8fj7+9P06ZN6dOnD/n5+TV+viOKKKWUn1JqilLq70qpfyml\nbEqpIdXs95RS6gul1Aml1Dml1AGlVIZSyrPPxwiH6dWrF9988w27du0iPz+f/Px87rnnHhISEsjP\nz/eIK8GuJTk5mZUrV7J+/XratGljdhxhh3p5Cu9ScjpPNCQ7d+6kZ8+etGnThilTplBWVsa8efPo\n0aMH27Zto1OnTjU6jgNO5zUHUoHDwC6gx1X2uxs4AKwE/gOEAM8ADyulwrXWJ+zsV9Qzfn5+V8wJ\n8vPzo1mzZpW3CJk4cSI//vgjixcvNiNirSUlJbFs2TJWrVqFn58fVqsVAH9//6ve09JsRUVFfPfd\nd5UjfwcOHCA/P5+AgABat25tcrorOTWv1tot2/33339BO1BxcbHu2bNnmWEY54AeZv/9pElzdAN0\nfHy8btasmf7Pf/6jKxw/flzfdNNN+tFHH9X2mjZtmgY0MMnOLN5Ay4tfRwE2YEgNnxt5cf8XzX5N\npbm2XXzvdV5enr6W2NhYPW7cuMrvExMTdWxs7DWfc7msrCxdk76cSSmlDcO4oi1evPiKfd0hr9Za\nf/HFF9XmHjZsWJX9PC2v1lrn5eVV/HsXqWvw81rvR6AqyEiUaAi++uor+vbty80331y5rVWrVnTv\n3p1PP/2U4uJiGjduXOPj1XYkSmtdCpy0L32lwxf/vPmae4kG6x//+EeV7zMzM01KUjc2m83sCHbr\n3r27R+V2Zt56PQfqcjInStR3v/zyC76+vldsb9y4MefPn2f37t12H9MVE8uVUgFKqRZKqXuATMp/\nC8xxRl9CCOEIDaqAAimiRP122223sWXLlorTIQCUlpaydetWAH788cdaHdcFRdSPgBXYBkQDY7TW\nUkAJIdxWgyugQIooUX8988wz7Nu3j+HDh1NQUMDu3bsZPHgwJ06Uz8Wuyz22nFxExQN9gfHAD4Cf\ng48vhBAO1WDmQF1O5kSJ+uipp57CarXyxhtvsHjxYpRS3HPPPbz44otMnz6dJk2a1On4zlpsU2u9\n4eKXa5RSq4DdSqmzWut5jji+8CwWi4WCggKn9lGxYKUr+nIEyet8Bw8etO8JNZlpbkZz9FV4VyNX\n50mrLw3QJSUlWmutT506pTdt2qR3796ttdZ64sSJ2jAMXVBQoB3BnqvzsPMqvIvP2QRsMvs1leba\nBkQbhlHxs+X05sq+JK9ntIuZo2vy89pgR6AqyEiUqE9SU1N55pln6NSpE926davc/vnnn/OrX/2K\n0NBQh/Rz+UgUsAwYDrQDDgEZWuv9dejCF7ihTiGFJzpvs9nIysqqXOPJWSwWC6mpqS7pyxEkr/MV\nFBSQkJAAcL5GTzD7N46rNVeNQFWQkShpnt4eAB3g5aW9DENnZmbqCsuXL9dKKZ2WlqYdrWIkSoHN\nH0q7Q5k/lCooAxK11nCVESjAC7j58r8H0AUoBTLNfk2lubZRw3WgHMFd1imqKcnrfLIOVC3JSJTw\ndBuAdWVlDAWGDx/O8ePHOXDgAIsWLeKhhx5izJgxDu9z0KBBTEpNZYTWajY08gVKwBgDLIAMpdRd\nlI8mAfRXSlUs/TuH8otYjiilsoH/A4qAzkAi5auST3N4YCGEcJAGeRXe1cjVecLTtQfuANCaSZMm\n8dVXXzFjxgw+/vhjDMPxH/eMjAxuNgzm8N8qyZfy6kiBAsYBIyn/rW4g8KeL7RagGHiP8hGqSRef\n9hCwBLhHa/29wwMLjzNlyhQMw6jSLr+1i6eaOXMmXbp0oWnTpgQGBjJw4ED27dtndqzrmjt3LiEh\nIfj6+hIdHc327dvNjlQtZ7++UkBdRooo4cnaA2uBHobBo48+yv/93//xwgsv0KiRcwabDx06RLjW\nXH7XLl/ggfLTeNlaa69q2g9a61Kt9Xit9d1a61u01j5a6/Za65Fa6x+cElh4pDvvvBOr1cqJEyc4\nceIEX375pdmRHCI3N5fRo0ezdetW1q1bR2lpKb17967TciPOlp2dTUpKClOmTGHnzp2Eh4fTp08f\nCgsLzY52BWe/vnIKrxpyOk94shJgl1KMbNfO6X21a9eOtUpRwn9HoCoy7Cz/8pDTQ4h6r1GjRrRo\n0cLsGA5nsViqfL9o0SJatmxJXl4eMTExJqW6trS0NEaOHMmQIUMASE9P57PPPiMjI4MXX3zR5HRV\nOfv1lRGoq5CRKOGJSoAxwM9aM2LECKf3N3z4cE5rzZiLfVdkGA2cKf/3ZaHTQ4h6b//+/QQHB9Oh\nQwcSEhI4cuRI5WPDhg0jNjbWxHSOc+rUKZRSBAQEmB2lWqWlpeTl5dGzZ8/KbUopevXqxebNm01M\nVjOOfn2lgLoGKaKEJ+kOBHt5kWkYLFi4kI4dOzq9z06dOrFg4UIyDYNfeXkRZxjlGZRCl8+BetLp\nIUS9Fh0dzaJFi1izZg3p6ekcPHiQ+++/n6KiIgCCgoJo27atySnrTmvN2LFjiYmJcds5XoWFhZSV\nlREYGFig1u+XAAAgAElEQVRle2BgYOXdDtyVM15fOYV3HXI6T3iKjcAfxo1j5MiRLimeKiQmJhIT\nE8PChQs5dOgQI9u1Y8SIEWRnZzt8xXLR8PTp06fy6zvvvJMuXbrQtm1bPvzwQ4YNG8aMGTNMTOc4\nSUlJfPvtt5UreAvHcsbrKwVUDUgRJTzF1KlT8fG5fEq383Xs2JGZM2dW2eas276Ihs3f35/bbruN\n7777zuwoDpOcnIzFYiE3N5egoCCz41xV8+bN8fLywmq1VtlutVpp1aqVSamuz1mvr5zCqyE5nSeE\n/Zx8A2LRAJ09e5bvv//erQsNeyQnJ7Ny5UrWr19PmzZtzI5zTd7e3kRFRZGTk1O5TWtNTk5OlTsf\nuBNnvr5SQNlBiigh7CdFlKiLF154gY0bN3L48GG++uorBg4cSKNGjXjiiScAmDhxIkOHDjU5Ze0k\nJSWxZMkSli5dip+fH1arFavVyrlz58yOdlXjx4/nvffe4/3332fPnj08++yzFBcXk5iYaHa0Kzj7\n9ZVTeHaS03lC2E9O54naOnr0KE8++ST/+te/aNGiBTExMWzZsoVmzZoBcPz48SpX5XmS9PR0lFL0\n6NGjyvbMzMzKZQLczaBBgygsLGTSpElYrVYiIiJYs2aNWy4z4ezXVwqoWpAiSgj7SRElamPZsmXX\nfDwzM9NFSRzPZrOZHaFWkpKSSEpKMjvGdTn79ZVTeLUkp/OEsJ+czhNC1BdSQNWBFFFC2E+KKCFE\nfSAFVB1JESWE/aSIEkJ4OpkD5QAyJ0oI+8mcKHE1FouFgoICp/ZRsaCiK/pyBMnrfAcPHrRrf6W1\ndlKUunnggQcubNy40cvsHPYoKSmhX79+tvXr15fabDYpooRLKaV0SUmJKQtp1sX06dN55ZVXACZL\nEdWwKaWiDcPY7KrJ1YZheNREbsnrfBczd9Vab7nevjIC5UAyEiWE/WQkSlzivM1mIysri7CwMKd2\nZLFYSE1NdUlfjiB5na+goICEhASA8zXZXwooB5MiSgj7SRElLhUWFkZkZKRT+6g4reSKvhxB8rof\nmUTuBDKxXAj7ycRyIYQnkQLKSaSIEsJ+UkSJ6hw7dozBgwfTvHlzGjduTHh4OF9//bXZsRwiNzeX\n/v37ExwcjGEYrFq1yuxIV5Wenk54eDj+/v74+/vTrVs3Vq9ebXas65o7dy4hISH4+voSHR3N9u3b\nHXJcKaCcSIooIewnRZS41KlTp7jvvvu48cYbWbNmDQUFBbz55pvccsstZkdziKKiIiIiIpg3bx5K\nKbPjXFPr1q2ZNWsWX3/9NXl5ecTFxTFgwAC3vsouOzublJQUpkyZws6dOwkPD6dPnz4UFhbW/eBa\na7ds999//wVdTxQXF+uePXuWGYZxDuih3eD1lVb/GqDXrl2rlVJXNMMw9NatW7Wr7N+/X//v//6v\n/tWvfqUbN26sQ0ND9Z/+9CddXFxc42NMmzZNAxqYpGv3evgBU4C/A/8CbMCQy/ZRQCKwEvgBOAt8\nA7wM3FibfqXV6Wc4EtB5eXmVPwcvvfSSfuCBB7SjZWVl6cv7MptSSq9cubLax9wxr9ZaBwQE6IyM\njCu2u0vee++9V48ZM6bye5vNpoODg/WsWbOu2DcvL6/i35xIXYOfVxmBcgEZiRKuNnbsWLKysirb\nBx98QMeOHV3S99GjR/nNb37Dtm3bGD16NLNnz6Zbt25MnjyZJ598ssbHccBIVHMgFQgFdlH+D+Pl\nGgMZF/edDzwPbKW88LLUok/hYJ988gn33HMPgwYNIjAwkMjISBYsWFBln1dffZWQkBCTEjZMNpuN\n5cuXU1xcTNeuXc2OU63S0lLy8vLo2bNn5TalFL169WLz5s11Pr5checicnWecKWYmBgeeeQRU/p+\n//33+fnnn9m8eTOhoaEAPPXUU5SVlfHBBx9w+vRp/P39a3SsOl6ddwxopbU+qZSKAqqb+HAe6Kar\nrvmyUCl1GHhVKRWntf6HHX0KBztw4ADz588nJSWFl19+mW3btjFmzBhuvPFGBg8eDECLFi1c9gtC\nQ7d79266du3KuXPnuOmmm1ixYkXl59zdFBYWUlZWRmBgYJXtgYGB7N27t87HlxEoF5KRKOFKZ8+e\npayszOX9njlzBoCWLVtW2d6qVSsMw+CGG26w63i1HYnSWpdqrU/WYJ/qFsxbQfnpPc9YwKYes9ls\nREVFMXXqVMLDw3n66ad5+umnSU9Pr9xn1KhRfP755yambDhCQ0PJz89n27ZtPPfccwwZMoQ9e/aY\nHcsUUkC5mBRRwhWGDRtG06ZN8fHxIS4ujry8PJf13aNHD7TWDB8+nPz8fI4ePUp2djbp6ek8//zz\n+Pr62n1MEyaWB1380wEzTUVdBAUFXbEQY1hYGD/88INJiRq2Ro0a0b59e+6++26mT59OeHg4s2fP\nNjtWtZo3b46XlxdWq7XKdqvVSqtWrep8fCmgTCBFlHAWb29vHn30UWbPns2qVauYPn06u3fv5oEH\nHiA/P98lGfr06cPUqVP5/PPPufvuu2nTpg1PPvkkY8aM4f/9v/9X6+O6uIh6EThN+QR0YaL77rvv\nitMte/fupW3btiYlEpey2Wz88ssvZseolre3N1FRUeTk5FRu01qTk5NDt27d6nx8mQNlEpkTJZwh\nOjqaHj16VH7/P//zP/zud7+jc+fOTJgwAYvFNfOi27VrR/fu3Xn00UcJCAjgs88+Y/r06bRq1Yqk\npKRaH9cVK5YrpSYCccBzWuufHX18YZ9x48Zx3333MXPmTAYNGsTWrVtZsGAB7733XuU+c+fOZcWK\nFaxbt87EpLVTVFTEd999V3EVIgcOHCA/P5+AgABat25tcrqqJk6cSN++fWnTpg1nzpxhyZIlbNiw\ngbVr15od7arGjx9PYmIiUVFRdOnShbS0NIqLi0lMTKz7wWtyqZ4ZrT4tY3AtssSBNEc1QJeUlOjq\nPPHEE9rHx0fbbLZqH3ekZcuW6caNG+tjx45V2T5s2DDdpEkT/e9//7vOfdi7xAEQRTXLGFSz3/8C\nZcBfanJcaQ7/Gb5iGQOttf7ss8/0XXfdpX19ffXtt9+uFy5cWOXxV199VYeEhGh7uMtl9l988UXl\nUiOXtmHDhlXZzx3yjhgxQoeEhGgfHx8dGBioH3zwQZ2Tk1Ptvu6Qt8LcuXN127ZttY+Pj46Ojtbb\nt2+vdj97lzGQESiTyUiUcIXWrVtz/vx5ioqKaNKkiVP7mj9/PpGRkQQFBVXZ3r9/fxYvXszOnTuJ\ni4urUx/OGIlSSj0ILAY+AZ6r6/GE4zz00EM89NBDV3188uTJTJ482YWJHKd79+7YbDazY9TI5ctH\neIqkpKQ6jXxfjcyBcgMyJ0o4SmpqKvv3779i+/fff4+Pj4/Tiqf9+/czYcIEnnjiCXbv3k1RUdEV\n+5SWlgJw4cIFh/R5+ZwopVQnpdRMpdSyi392qumxlFL3An8DtgH/q7X2jP/RhBCmkQLKTUgRJerq\nAWDBW28RFhrKokWLKrfn5+fzySef0KdPH6f0m5mZSVhoKO++8QbWDz/kzL//zc6dO5k1a1aV/ZYu\nXYphGHTu3NlhfVcpomCvP/yhOwzyhz8o2KOUSrzeMZRSYcCnwAGgn9baPWfECiHcipzCcyNyOk/U\nxQYg1mbje2D48OEcO3aMn376iffee48mTZowc+ZMh/e5f/9+nn7qKYbbbMwGfIF1QB/gj3/8IydO\nnOC2227jk08+Yc2aNTz99NMOuXz4UoMGDWJSaiojtFazoZEvUALGGGABZCil7roYDaC/UqpiZu4c\nyuc7rAFuBl4H/uey+5F9r6tfJ0oI0cBJAeVmpIgSdfE74APgqNZMmjSJVq1a8eijjzJp0iTat2/v\n8P4yMjLwV4o5gM/Fbb0oL+biKJ8zUVpaSkhICDNmzOCFF15wSoabDYM5ZWWVGXwpr44WgtIwjvJC\nSQMDLzYof6kUEHzx+9eqOfxiQAooF7NYLE6/Qe2mTZtc1pcjSF7nO3jwoH1PqMlMczNaQ7kK72rk\n6jxpdjfQFS3WMPTjjz+une3xxx/XsYahL+3bnTJ0hwvAMm32eyOtRg2INgyjoth1enNlX5LXM9rF\nzNE1+XmVESg3JSNRorZKgF1KMbJdO6f31a5dO9YqRQn/PUfmThl2ln95yOkhhKOct9lsZGVlXbH6\nuKNZLBZSU1Nd0pcjSF7nKygoICEhAcrvkXl9Zv/GcbXW0EegKshIlLQaN9DFoJ8C7WUYev/+/drZ\n9u3bp70MQz91sW93yjACtCpf+6mjNvu9kVajxlXWgXIGd1qnqCYkr/PZuw6UXIXn5uTqPFFT3YFg\nLy8yDYMFCxe65O70nTp1YsHChWQaBr/y8iLOMNwng1Lo8jlOTzo9hBCiwZECygNIESVqYiMwYtw4\n9uzd65jbFNRQYmIie/bu5ZkXXiBw0CBGvvCCW2TYu2+fq29ALJwgJCQEwzCuaKNHjzY7mkPk5ubS\nv39/goODMQyDVatWmR3pqmbOnEmXLl1o2rQpgYGBDBw4kH379pkd65qc+frKHCgPIXOiRE1MnToV\nHx+f6+/oYB07dnTKMgl1zeCKe+cJ59qxYwdlZWWV33/zzTf07t2bQYMGmZjKcYqKioiIiGDEiBE8\n8sgjZse5ptzcXEaPHs0999zDhQsXmDBhAr1796agoABfX9/rH8AEznx9pYDyIFJECWE/KaI8W7Nm\nzap8/8knn9ChQwfuv/9+kxI5Vnx8PPHx8QAV88Dc1uU3I1+0aBEtW7YkLy+PmJgYk1JdmzNfXzmF\n52HkdJ4Q9rv8ti9m5xG1U1paypIlSxgxYkTltmHDhhEbG2tiqobr1KlTKKUICAgwO4oppIDyQFJE\nCWE/KaI834oVKzh9+jRDhw6t3BYUFETbtm1NTNUwaa0ZO3YsMTEx3H777WbHMYWcwvNQcjpPCPvJ\n6TzPlpGRQd++favcDmjGjBkmJmq4kpKS+PbbbytXHG+IZATKg8lIlBD2k5Eoz/TDDz+wbt06nn76\nabOjNHjJyclYLBa++OILgoKCzI5jGimgPJwUUULYT4ooz5ORkUFgYCAPPfSQ2VEatOTkZFauXMn6\n9etp06aN2XFMJQVUPSBFlBD2kyLKc2itWbRoEYmJiRhG1f+2Jk6cWGVOlKcpKioiPz+fXbt2AXDg\nwAHy8/M5cuSIycmulJSUxJIlS1i6dCl+fn5YrVasVivnzp0zO9pVOfP1lQKqnpAiSgj7SRHlGdat\nW8eRI0cYNmzYFY8dP37cLYuNmtqxYwd33303UVFRKKVISUkhMjKSyZMnmx3tCunp6fz888/06NGD\nW2+9tbJ9+OGHZke7Kme+vjKJvB6RieVC2E8mlru/Bx98sMpimpfKzMx0cRrH6t69OzabzewYNeIp\nOS/lzNdXRqDqGRmJEsJ+MhIlhLCXFFD1kBRRQthPiighhD2kgKqnpIgSwn5SRAkhakrmQNVjMidK\nCPvJnCjzWSwWCgoKnNpHxQKQrujLESSv8x08eNCu/ZW73rzwgQceuLBx40Yvs3PUByUlJfTr18+2\nfv36UpvNJkVUPaWU0iUlJfj4+JgdpV6YPn06r7zyCsBkKaJcQykVbRjGZldNVjYMw6MmRkte57uY\nuavWesv19pURqAZARqKEsJ+MRJnivM1mIysri7CwMKd2ZLFYSE1NdUlfjiB5na+goICEhASA8zV6\ngtbaLdv9999/QQuHKi4u1j179iwzDOMc0EO7wfsszXEN0F9++aUeNWqUvuOOO7Sfn59u06aNHjRo\nkN63b58207Rp07RSSt91110u6S8xMVErpapthmHoY8eO1fhY06ZN04AGJmn73xM/YArwd+BfgA0Y\nUs1+vwHmATsu/uNdZm9f9aEBkYDOy8ur8ftTW1lZWdpVfTmC5HW+vLy8is96pK7Bz6uMQDUgMhJV\n/7355pts3bqVxx57jM6dO3PixAn+/Oc/ExkZydatW025a/qPP/7IzJkzadKkicv6fPbZZ3nwwQer\nbNNaM3LkSNq3b2/X/bvqOBLVHEgFDgO7gB5X2e8hYDjwT+B74DY7+hBCmECuwmtg5Oq8+u3555/n\n8OHDvP322wwfPpyJEyeSm5vLhQsXeO2110zJlJKSQteuXYmKinJZn/feey9PPvlkldauXTuKi4v5\n/e9/b/fx6nB13jGgldY6BHgRUFfZbx7gr7XuAqyzO2A9ZrPZSE1NpX379jRu3JiOHTtWvBf1Qm5u\nLv379yc4OBjDMFi1apXZka7Jk/LOnDmTLl260LRpUwIDAxk4cCD79u1z2PGlgGqApIiqv+69914a\nNao6sNyxY0fuuOMOU66E2bhxI3/72994++23Xd735ZYsWYJhGDzxxBO1en5tiiitdanW+mQN9vtJ\na/1LrYLVc6+99hp/+ctfmDdvHnv27OH111/n9ddf55133jE7mkMUFRURERHBvHnzUOpq9bX78KS8\nubm5jB49mq1bt7Ju3TpKS0vp3bs3JSUlDjm+nMJroOR0XsNitVq58847XdqnzWZjzJgxPP3009xx\nxx0u7ftyFy5c4KOPPuK+++6r0x3kZWK5623evJkBAwYQHx8PQJs2bVi6dCnbtm0zOZljxMfHV/7d\ntHbPq+Iv5Ul5LRZLle8XLVpEy5YtycvLIyYmps7HlxGoBkxGohqGrKwsfvzxRx5//HGX9jt//nx+\n+OEHpk6d6tJ+q7N69Wr+9a9/1er03eVksU3X6tatGzk5Oezfvx+A/Px8Nm3axEMPPVS5z6uvvkpI\nSIhZEYWHOHXqFEopAgICHHI8GYFq4GQkqn7bs2cPycnJ3HfffQwZMsRl/f773/9m8uTJTJo0yWH/\nWNXF0qVLueGGG3jssccccjwZiXKdP/7xj/z888+Ehobi5eWFzWZj+vTpVX4haNGiBR07djQxpXB3\nWmvGjh1LTEyMwy6mkREoISNR9ZTVauXhhx/mlltu4aOPPnLpfIWXX36ZZs2akZyc7LI+r6aoqIhV\nq1YRHx/PLbfc4rDjykiUa2RnZ7N06VKWL1/Ozp07Wbx4MW+88QYffPBB5T6jRo3i888/NzGlcHdJ\nSUl8++23LF++3GHHlBEoAchIVH2RmprKM888Q2BgIPHx8fz88898+eWXtGrVyqn97t+/n4yMDA4d\nOoS/vz8LFixg9uzZ/Pjjj0D5b3/nzp2jtLSUw4cP07RpU4cWM5dnaNeuHcOHD6dTp06sWLGCkpIS\nh5y+u9zlI1HAMsqXI2gHHAIytNb7Hd5xA/Liiy8yceLEytHDO+64g0OHDjFz5kwGDx5scjrhCZKT\nk7FYLOTm5tq1hMn1SAElKkkR5dkeADLS0njrzTfp2KkTx44dIycnh1//+tdO7TczM5Onn3oKf6UI\n15rPgDKbjdGjRzN69Ogr9m/fvj3PP/88b731ltMyrFWKN15/nQULF5KdnU2TJk3o16+fw/q7VJUi\nCl5tCmURYOwC28/wolJqhNZ6kVM6bwCKi4uvGD31xFuECHMkJyezcuVKNmzYUKcLSKojp/BEFXI6\nz3NtAH4oK+NXWrNv3z7mzJlDly5dnNrn/v37efqppxhus3G0rIx/2Gx8a7PxIOULHs2fP5+PP/6Y\njz/+mDvuuIO2bdvy8ccfM2LECKdmOFpWxjCbjRHDh5OTk8Mjjzzi1HsEDho0CEMpngJ1HBp9AcZx\naDQCDAULlVIyQaeW+vXrx/Tp07FYLBw+fJgVK1aQlpbGI488UrnP3Llz6dWrl4kpa6+oqIj8/Hx2\n7doFwIEDB8jPz+fIkSMmJ6ueJ+VNSkpiyZIlLF26FD8/P6xWK1arlXPnzjnk+DICJa4gI1Ge62Xg\nCOANrFixghtuuKHK444+jZWRkYG/UswBKsqTXwErgWDD4NChQ4wcORKAtLQ0lFIOHwmqLoMvMAdY\nApy7cMEpp+8uz3CzYTCnrOyKDFnAOVgI/N/Fh/orpVpf/HqO1vqMUqoNUHE+6h4ApdTLF78/rLXO\ncupfwI298847pKamMmrUKE6ePMmtt97Kc889R2pqauU+hYWFHDhwwMSUtbdjxw5iY2NRSqGUIiUl\nBYChQ4eSkZFhcroreVLe9PR0lFL06NGjyvbMzEyHXFQjBZSolhRRnimf8pGfUuDTTz/ls88+q/K4\nowuJQ4cOEa41l4/t+AIRWnPo0KEq250xkf1aGby1Rvv40LNnT4f3W9MMlI/03w/EUH6frYEXG8AH\nwBkgBJh68fEKFVf2baC8DmuQ/Pz8eOutt655ynfy5MlMnjzZhakcp3v37h51OtKT8jo7p5zCE1cl\np/M8z3rgLHCLlxcvvfQSZWVlVZqjtWvXjnyluHxd3xJgl1K0a9fuv9nWryc/P9+lGby8vBg7dqzT\nr0C8VoYboAyYpbX2qqb9AKC13qC1Nq6yT5xTwwshakUKKHFNUkR5lhJgDPCz1g6dZ3Q1w4cP57TW\njLnYt2SommE0cKb839mFTg8hhHApKaDEdUkR5Rm6A8FeXmQaBgsWLnTJwoKdOnViwcKFZBoGv/Ly\nIs4wJMOlGZRCl59VfdLpIYQQLiVzoESNyJwo97cR+MO4cYwcOdKlqzInJiYSExPDwoULyyeNt2vH\niBEjJMPFDNnZ2bJiuZ0sFovTb369adMml/XlCJLX+Q4ePGjX/spdbwb4wAMPXNi4caOX2TlEVSUl\nJfTr18+2fv36UpvNJkWUG1FK6ZKSEqderi9qZ/r06bzyyisAk6WIujqlVLRhGJtdNUnZ09aTkrzO\ndzFzV631luvtKyNQwi4yEiWE/eTeeTV23mazkZWVRVhYmFM7slgspKamuqQvR5C8zldQUEBCQgLA\n+ZrsLwWUsJsUUULYT4qomgsLCyMyMtKpfVScVnJFX44ged2PTCIXtSITy4Wwn9yAWIj6QwooUWtS\nRAlhPymi7Hf27FnGjh1Lu3btaNy4MTExMezYscPsWA6Rm5tL//79CQ4OxjAMVq1aZXak65o7dy4h\nISH4+voSHR3N9u3bzY5UrfT0dMLDw/H398ff359u3bqxevVqhx1fCihRJ1JECWE/KaLsM2LECHJy\ncliyZAm7d+/mwQcfpFevXhw/ftzsaHVWVFREREQE8+bNc/qCr46QnZ1NSkoKU6ZMYefOnYSHh9On\nTx8KCwvNjnaF1q1bM2vWLL7++mvy8vKIi4tjwIABDrsqUAooUWdSRAlhPymiaubcuXP87W9/4403\n3uC+++6jffv2TJ48mY4dOzJ//nyz49VZfHw8f/rTnxgwYADuelX8pdLS0hg5ciRDhgwhNDSU9PR0\nGjdu7Hb3wQN4+OGHiY+Pp0OHDnTs2JFp06bRpEkTtmy57gV2NSIFlHAIKaKEsJ8UUdd34cIFysrK\nuPHGG6ts9/X15csvvwTg1VdfJSQkxIx4DUppaSl5eXlV7i2plKJXr15s3rzZxGTXZ7PZWL58OcXF\nxXTt2tUhx5QCSjiMFFFC2E+KqGtr0qQJXbt2ZerUqRw/fpyKZQ42b95ceQqvRYsWLl00taEqLCyk\nrKyMwMDAKtsDAwM5ceKESamubffu3dx0003ceOONJCUlsWLFCkJDQx1ybCmghENJESWE/aSIuras\nrCy01gQHB+Pj48M777zDk08+iWGU/xc2atQoPv/8c5NTCncUGhpKfn4+27Zt47nnnmPIkCHs2bPH\nIceWAko4nBRRQthPiqirCwkJYf369RQVFXHkyBG2bNnC+fPnad++vdnRGpTmzZvj5eWF1Wqtst1q\ntdKqVSuTUl1bo0aNaN++PXfffTfTp08nPDyc2bNnO+TYUkAJp5AiSgj7SRF1bb6+vgQGBvKf//yH\nNWvW8Nvf/tbsSA2Kt7c3UVFR5OTkVG7TWpOTk0O3bt1MTFZzNpuNX375xSHHkgJKOI0UUULYT4qo\nK61du5Y1a9Zw6NAhPv/8c+Li4rj99ttJTEwEytcl6tWrl7kha6moqIj8/Hx27doFwIEDB8jPz+fI\nkSMmJ6ve+PHjee+993j//ffZs2cPzz77LMXFxZXvhTuZOHEiubm5HD58mN27dzNhwgQ2bNhQcbuW\nOpNbuQinktu+CGE/ue1LVadPn2bChAn8+OOPBAQE8OijjzJt2jS8vMrvN19YWMiBAwdMTlk7O3bs\nIDY2FqUUSilSUlIAGDp0qFsuDTBo0CAKCwuZNGkSVquViIgI1qxZQ4sWLcyOdoWTJ08ydOhQjh8/\njr+/P507d2bt2rXExcU55PhSQAmnkyJKCPtJEfVfjz32GI899thVH588eTKTJ092YSLH6d69Ozab\nzewYdklKSiIpKcnsGNe1YMECpx5fTuEJl5DTeULYT07nCeG+3L6AOn/+PC+99BLBwcE0btyY6Oho\n1q1bZ3YsUQuOKKKUUn5KqSlKqb8rpf6llLIppYY4Ia7HcofPTFFREZMnT6Zv3740a9YMwzB4//33\nXdb/t99+y6BBg+jQoQN+fn60aNGC7t278+mnn7osQ4Wvv/6a/v3706xZM/z8/Ljrrrt45513avz8\nuhZR9nxmlFKhSqnVSqkzF/d9XynV3N4+hWgI3L6AGjp0KG+//TaDBw9mzpw5NGrUiIceeoivvvrK\n7GiiFhxQRDUHUoFQYBfg/vc+cDF3+MwUFhYydepU9uzZQ0REhMvv8XX48GHOnj1LYmIic+bMYdKk\nSSil6N+/v9OH9S+1du1aunXrVjlnZM6cOfTr14+jR4/adZw6FlE1+swopYKBXKA98EfgDeBhYK1S\nSqZ7CHEZt/5QbNu2jezsbN58803GjRsHwODBg7nzzjt58cUXK5fxF56ljnOijgGttNYnlVJRgHve\nBtwk27dvd4vPzK233sqJEydo2bIleXl5/OY3v3FJvxX69u1L3759q2xLTk4mMjKSt956i6eeesrp\nGc6cOcPQoUPp168fH330UZ2PV4c5UTX9zLwM+AIRWusfAZRS24HPgUTAZZWnxWJx2A1fr2bTpk0u\n67C2G20AABw8SURBVMsRJK/zHTx40L4naK3dst1///0XXnjhBe3t7a3PnDmjLzVz5kxtGIY+evSo\nFp6ruLhY9+zZs8wwjHNAD23nzwgQBdiAIfY+tz42QI8bN87tPjM7duzQSim9ePFil/d9uX79+umg\noCCX9DV//nxtGIbeu3ev1lrroqIibbPZ6nzcadOmacpHkSZpB35mgBPA8mq27wHW2ttXbRoQbRhG\nxd/P6c2VfUlez2gXM0fX5OfVrUegdu3axW233UaTJk2qbO/SpUvl48HBwWZEEw4gV+c53j//+U/5\nzFyiuLiYkpISTp8+zcqVK/n73//OE0884ZK+c3JyaNq0KUeOHKF///7s27cPPz8/Bg8eTFpa2hU3\nx60pZ1ydp5S6FWgJ7Kjm4W1A32q2O8P5invdhYWFObUji8VCamqqS/pyBMnrfAUFBRVrRJ2vyf5u\nXUAdP36coKCgK7YHBQWhtebYsWMmpBKOJEWUY504cUI+M5dISUnhL3/5CwCGYfC73/2OP//5zy7p\ne//+/ZSWljJgwACefvppXnvtNb744gvmzJnD6dOnWbJkSa2P7YQiquKH5ng1jx0HApRS3lrr0jr2\nUyNhYWFERkY6tY+K00qu6MsRJK/7cetJ5CUlJdX+lubj41P5uPB8ssSB48hnpqpx48axbt063n//\nfR566CHKysocdhuH6zl79iwlJSUkJiaSlpbGb3/7W95++21GjhzJ8uXL+f777+t0fAcvceB78c/q\nXpxzl+0jhMDNCyhfX99q/7E7d+5c5eOifpAiyjHkM1PVbbfdRlxcHAkJCaxatYozZ87Qv39/l/Rd\n8Vo//vjjVbY/+eSTaK3ZvHlznftwYBFVUVlXd17R57J9XC43N5f+/fsTHByMYRisWrXqin0mTZrE\nrbfeSuPGjXnwwQf57rvvTEhaO3PnziUkJARfX1+io6PZvt0zro157bXXMAyD8ePHmx3lmpz1+rpt\nAVVQUGCcPXuWw4cPX/HY8ePlo8y33nqrq2MJJ6quiFJKdVJKzVRKLTM7nycw6zOzf/9+JkyYwBNP\nPMGECRPYv3+/U/qpa4ZHH32U7du3Oy3fpRmKiooACAwMrLJPy5YtAfjPf/7jkD4vL6Jq+ZmpOHV3\n5fnf8m3/dtXpu+oUFRURERHBvHnzql0SY9asWbzzzju8++67bNu2DT8/P/r06cP58zWaymKq7Oxs\nUlJSmDJlCjt37iQ8PJw+ffpQWFhodrRr2r59O++++y7h4eFmR7kmp76+rriyolZXY7RDK2+lAT1/\n/vwqV6FMnz5drsKrxyquzlNKlQJl3EgpbSnTdlxR1BDbA6B9lOs/MxkZGdrLMHSAl5f+/+3dfVBU\n9+Hv8c/ZdW2IMRjEAiFo9PZB0hkhkCLSpAalkUxucDq9YdJeBojpHTubh1romI4zYuhNbMwfMjqi\ntkmU9Gqid+I4oZUEG7BqHJ/YwBqnhNbBm6LVjdtprIFGkT33D2X7IzyeurvnYN6vGUY5e+Z8P3t2\nVz9zzvfsyb/+p9vlMrdu3WqaZmyuwhstQ79169aZLpfLPH78eNQz9L8Wzz333ID1mpubTcMwzDff\nfDOi4/dfnWdIoXipd76sfWYkBTT8VXh/+OLyaPxIypJk+ny+YZ+nYRjm22+/PWBZSkqKuXbt2vDv\nFy9eNG+55RZz586dw25n27Zt5mhjxcLcuXPNZ599Nvx7KBQyU1NTzTVr1gxYzyl5TdM0L126ZH7j\nG98wm5qazAcffND82c9+Nmgdp+Qd6/41TdP0+Xz9V+NlmWN4vzr2CJTKJfN/mpKu3Xen/3DslStX\nVFdXp9zc3C/V1URfJnFxcaqpqXGZMicoSy79XBP0hIPfqw6xX1KjGdvPzF/+8hf9rx//WEtCIZ3p\n61Pz9T+fCIX04yefjMlplKEytA6R4erVq3r99dcVFxene+65J+oZ9l1/LV5es2bAfnjllVfk8Xj0\n4IMPRjRDcXGxXIahH0vGOWnCH62fYdgl6b9f/0JNSZJhGAslfUPS/41g1Ig6ffq0zp8/r4ULF4aX\n3X777Zo7d+6A06T5+flasmSJHRGH1dvbK5/PNyC7YRgqKCiIyCneaHnqqaf06KOPRuymvNES7f3r\n6KvwdLekdMlsN1VSUqInn3xSdXV1+vjjj7V161a70yGK3njjDblvdavv4T7JM/AxwzCekjRFUv8/\n9EWGYaRd//t60zQvxS6ps3xX0g8k7TJj85nZsmWL4g1D6/XviTJxktZL2m6aevLJJ/Wtb31LklRf\nX6+uri5J0rPPPqvJkydHLcMySZ9K8pimnn32WX3nO9/R9u3b1dHRobVr1+rWW2+NyNgjZciVVCbp\ndUlFRUV65plntG/fPu3atUsrVqxQcnJyxDNMcbm0vq8vnKHfGD8zqyX9D0l/NAxjnaTJkn4uyS+p\nLqJhI+j8+fMyDGPQqdKkpCSdP38+/PuMGTOGvELVTsFgUH19fUNm7+josCnVyHbs2KG2tja1tAz1\njRfOEu396+wCJV3732C9dPToUR09elRut9uMi4sLPfroo6bd0RA9n332mbvvzj7ji+Xpup9Lmn79\n76ak71//kaT/I+lLW6Ak6Q1J/02x+cx89tln7nl9fcYX/8OOkxQyTR04cEAHDhyQJO3atUu7du2S\nJL388stX3W531DI8Lum16xneeecdvfvuu+H9UF1dbVZXV0dk7JEy6HqGP+raJd1er1cul0tf+cpX\nQrW1taHa2tqYZLhu1M+MaZpnDMOYL2mtpF/p2nfh/F7Sz00b5z9FSl1dnd0Rxr0zZ85o2bJleu+9\n9+TxDP2P85eJYwuUucr890zB/21jEDiOaZoz7c7gSOa1z8xESV02R5H+fe27HYqv/9jNLen/2R1C\nY//MmKbZrth9aWZEJCcnyzRNBQKBAUcaAoGA7r33XhuTjS4xMVFut1uBQGDA8kAgEPEjlJHg8/l0\n4cIFZWVl9c9ZU19fnw4cOKANGzbo8uXLMb/v5UiivX+ZVwIAGLdmzpyp5ORkNTU1hZf985//1NGj\nR5WXl2djstF5PB5lZ2cPyG6appqamhyZvaCgQB9++KHa2trk9/vl9/t13333qaSkRH6/31HlSYr+\n/nXsESgAAKRrX2Nw6tSp8FGPzs5O+f1+JSQkKC0tTcuWLdMLL7ygr33ta7r77ru1cuVK3XXXXVq8\neHF4G2VlZUpNTdXq1avtehpDqqioUHl5ubKzs5WTk6Oamhr19PSovLzc7miDTJo0adAFGJMmTdLU\nqVMde7uWaO5fChQAwNFaWlqUn58vwzBkGIYqKyslXStFW7Zs0fLly9XT06OlS5fq008/1QMPPKB3\n3nlHEydODG+jq6tLkZp3F0nFxcUKBoOqqqpSIBBQZmamGhsbNW3aNLujjYnTjjp9UTT3LwUKAOBo\n8+fPVygUGnGd559/Xs8///ywjzc3N0c4VeR4vV55vV67Y/xHnLxf+0Vr/zIHCgAAwCIKFAAAgEUU\nKAAAAIuYAwUAcJSGhga1t7dHdYxDhw7FbKxIIG/0nT592tL6Rv9loQAA2MkwjFyXy3V4tAnjkeJy\nuUadnO4k5I2+65nnmaZ5ZLR1OQIFAHCKK6FQSNu2bYv69wo1NDRo5cqVMRkrEsgbfe3t7SopKZGu\n3cZoVBQoAICjpKenKysrK6pj9J9WisVYkUBe52ESOQAAgEUUKACAox08eFBFRUVKTU2Vy+VSfX39\ngMd3796tRYsWKTExUS6XSydOnLAp6X+mtrZWM2fOVFxcnHJzc3X8+HG7Iw1p8+bNysjIUHx8vOLj\n45WXl6d3333X7lgjGu29cyMoUAAAR+vu7lZmZqY2btw45K1Duru79cADD+jll192/K1Fvmjnzp2q\nrKxUdXW1WltblZGRoUWLFikYDNodbZC0tDStWbNGH3zwgXw+nxYsWKDFixc7+iq70d47N4I5UAAA\nRyssLFRhYaEkaagrx69P/NXHH3885ONOVlNTo6VLl6q0tFTStaM8e/bsCd/jz0keeeSRAb+/8MIL\n2rRpk44cOeLYieKjvXduBEegAAA3vSeeeEL5+fl2xxigt7dXPp9PCxcuDC8zDEMFBQU6fPiwjclG\nFwqFtGPHDvX09GjevHl2x7EFR6AAADe9lJQUxx2dCgaD6uvrU1JS0oDlSUlJ6ujosCnVyE6ePKl5\n8+bp888/1+TJk7V7927Nnj3b7li2oEABAG56q1evtjvCTWH27Nny+/26ePGi3nrrLZWWlurAgQNf\nyhJFgQIAwAaJiYlyu90KBAIDlgcCASUnJ9uUamQTJkzQrFmzJEn33nuvjh07pnXr1mnTpk02J4s9\n5kABAG4a4+kqPI/Ho+zsbDU1NYWXmaappqYm5eXl2Zhs7EKhkC5fvmx3DFtwBAoA4Gjd3d06depU\neA5TZ2en/H6/EhISlJaWpn/84x/661//qrNnz8o0TX300UcyTVPJycnh+UUrVqzQ2bNn9frrr9v5\nVAapqKhQeXm5srOzlZOTo5qaGvX09Ki8vNzuaIOsWLFCDz/8sKZPn65Lly5p+/bt2r9/v/bu3Wt3\ntGGN9t65ERQoAICjtbS0KD8/X4ZhyDAMVVZWSpLKysq0ZcsW1dfX64knngg//sMf/lCStGrVKlVV\nVUmSzp07p66uLtuew3CKi4sVDAZVVVWlQCCgzMxMNTY2atq0aXZHG+STTz5RWVmZzp07p/j4eM2Z\nM0d79+7VggUL7I42rNHeOzeCAgUAcLT58+crFAoN+3hZWZnKyspG3MbWrVsjHStivF6vvF6v3TFG\n9eqrr9odwbLR3js3gjlQAAAAFlGgAAAALKJAAQAAWMQcKACAozQ0NET9BrWHDh2K2ViRQN7oO336\ntKX1Dad9tT0A4MvJMIxcl8t1OFqTfr/I5XJFbYJxNJA3+q5nnmea5pHR1uUIFADAKa6EQiFt27ZN\n6enpUR2ooaFBK1eujMlYkUDe6Gtvb1dJSYkkXRnL+hQoAICjpKenKysrK6pj9J9WisVYkUBe52ES\nOQAAgEUUKACAox08eFBFRUVKTU2Vy+VSfX19+LGrV6/queee05w5c3TbbbcpNTU1/G3Z40Vtba1m\nzpypuLg45ebm6vjx43ZHGtKvfvUr5eTk6Pbbb1dSUpK+//3v689//rPdsYYV7bwUKACAo3V3dysz\nM1MbN24cdLPgnp4etbW1adWqVWptbdXu3bvV0dGhxYsX25TWmp07d6qyslLV1dVqbW1VRkaGFi1a\npGAwaHe0QQ4ePKhnnnlGR48e1Xvvvafe3l499NBD+te//mV3tCFFOy9zoAAAjlZYWKjCwkJJ0hev\nHL/99tvV2Ng4YNmGDRs0d+5cnTlzRnfddVfMcv4nampqtHTpUpWWlkqSNm/erD179mjLli1avny5\nzekGamhoGPB7XV2dvvrVr8rn8+n++++3KdXwop2XI1AAgJvKp59+KsMwNGXKlPCy/Px8LVmyxMZU\ng/X29srn82nhwoXhZYZhqKCgQIcPH7Yx2dj07+eEhAS7o4xJpPNSoAAAN43Lly/rF7/4hX70ox/p\ntttuCy+fMWOGUlJSbEw2WDAYVF9fn5KSkgYsT0pK0vnz521KNTamaWrZsmW6//77dc8999gdZ1TR\nyMspPADATeHq1at67LHHZBiGNm7cOOCxuro6e0LdpLxer/70pz+Fv3Hc6aKRlwIFABj3+stTV1eX\nmpubBxx9cqrExES53W4FAoEBywOBgJKTk21KNbqnn35aDQ0NOnjwoOOO6g0lWnk5hQcAGNf6y1Nn\nZ6eampp0xx132B1pTDwej7Kzs9XU1BReZpqmmpqalJeXZ2Oy4T399NN6++23tW/fPk2fPt3uOKOK\nZl6OQAEAHK27u1unTp0KX4HX2dkpv9+vhIQEpaSk6Ac/+IHa2tr0+9//Xr29veEjOgkJCfJ4PJKk\nsrIypaamavXq1bY9j6FUVFSovLxc2dnZysnJUU1NjXp6elReXm53tEG8Xq/efPNN1dfXa9KkSeH9\nHB8fr1tuucXmdINFOy8FCgDgaC0tLcrPz5dhGDIMQ5WVlZKulaJVq1bpd7/7nQzDUGZmpqRrR3EM\nw9C+ffv03e9+V5LU1dUlt9tt23MYTnFxsYLBoKqqqhQIBJSZmanGxkZNmzbN7miDbN68WYZh6MEH\nHxywfOvWreGvYXCSaOelQAEAHG3+/PkKhULDPj7SY/2am5sjGSmivF6vvF6v3TFGNZb97CTRzssc\nKAAAAIsoUAAAABZRoAAAACxiDhQAwFEaGhrU3t4e1TH6v1AxFmNFAnmj7/Tp05bWN754Y0YAAOxg\nGEauy+U6HKvJyi6Xa1xNjCZv9F3PPM80zSOjrcsRKACAU1wJhULatm2b0tPTozpQQ0ODVq5cGZOx\nIoG80dfe3q6SkhJJujKW9SlQAABHSU9PV1ZWVlTH6D+tFIuxIoG8zsMkcgAAAIsoUAAARzt48KCK\nioqUmpoql8ul+vr6AY9XV1crPT1dt912mxISEvS9731Px44dsymtdbW1tZo5c6bi4uKUm5ur48eP\n2x1pWKO9Fk6yefNmZWRkKD4+XvHx8crLy9O7774bse1ToAAAjtbd3a3MzExt3LhRhmEMevyb3/ym\namtrdfLkSR06dEh33323HnroIf3973+3Ia01O3fuVGVlpaqrq9Xa2qqMjAwtWrRIwWDQ7mhDGu21\ncJK0tDStWbNGH3zwgXw+nxYsWKDFixdH7KpA5kABABytsLBQhYWFkqShrhx//PHHB/y+du1avfba\nazpx4oTy8/NjkvE/VVNTo6VLl4bvzbZ582bt2bNHW7Zs0fLly21ON9hor4WTPPLIIwN+f+GFF7Rp\n0yYdOXIkIhPbOQIFALhp9Pb26te//rWmTJmijIyM8PL8/HwtWbLExmSD9fb2yufzaeHCheFlhmGo\noKBAhw8ftjHZzScUCmnHjh3q6enRvHnzIrJNjkABAMa9PXv26PHHH1dPT4/uvPNO/eEPf1BCQkL4\n8RkzZiglJcXGhIMFg0H19fUpKSlpwPKkpCR1dHTYlOrmcvLkSc2bN0+ff/65Jk+erN27d2v27NkR\n2TYFCgAw7i1YsEB+v1/BYFCvvPKKHnvsMR07dkyJiYmSpLq6OnsDwhazZ8+W3+/XxYsX9dZbb6m0\ntFQHDhyISIniFB4AYNyLi4vTrFmzlJOTo1deeUUTJkzQa6+9ZnesESUmJsrtdisQCAxYHggElJyc\nbFOqm8uECRM0a9Ys3XvvvXrxxReVkZGhdevWRWTbFCgAwE0nFArp8uXLdscYkcfjUXZ2tpqamsLL\nTNNUU1OT8vLybEx284rk+4JTeAAAR+vu7tapU6fCV311dnbK7/crISFBU6dO1YsvvqiioiKlpKQo\nGAxqw4YN+tvf/qbHHnssvI2ysjKlpqZq9erVdj2NIVVUVKi8vFzZ2dnKyclRTU2Nenp6VF5ebne0\nIY30WqSlpdmcbqAVK1bo4Ycf1vTp03Xp0iVt375d+/fv1969eyOyfQoUAMDRWlpalJ+fL8MwZBiG\nKisrJV0rRZs2bdJHH32k3/72twoGg5o6daq+/e1v6/333x9wqXpXV5fcbrddT2FYxcXFCgaDqqqq\nUiAQUGZmphobGzVt2jS7ow1ppNdiy5YtNqcb6JNPPlFZWZnOnTun+Ph4zZkzR3v37tWCBQsisn0K\nFADA0ebPn69QKDTs47t27Rp1G83NzZGMFFFer1der9fuGGMy2mvhJK+++mpUt88cKAAAAIsoUAAA\nABZRoAAAACxiDhQAwFEaGhoidsPX4Rw6dChmY0UCeaPv9OnTltY3nH4zQADAl4NhGLkul+twrCYp\nu1yucTMhWiJvLFzPPM80zSOjrcsRKACAU1wJhULatm3bgK8giIaGhgatXLkyJmNFAnmjr729XSUl\nJZJ0ZSzrU6AAAI6Snp6urKysqI7Rf1opFmNFAnmdh0nkAAAAFlGgAACOdvDgQRUVFSk1NVUul0v1\n9fXDrvuTn/xELpdL69evj2HCG1NbW6uZM2cqLi5Oubm5On78uN2RRjTe8vZ76aWX5HK5VFFREZHt\nUaAAAI7W3d2tzMxMbdy4UYZhDLve7t27dfToUaWmpsYw3Y3ZuXOnKisrVV1drdbWVmVkZGjRokUK\nBoN2RxvSeMvb7/jx4/rNb36jjIyMiG2TAgUAcLTCwkL98pe/1OLFizXcleNnz57VT3/6U73xxhua\nMGH8TO+tqanR0qVLVVpaqtmzZ2vz5s269dZbHXdfuX7jLa8kffbZZyopKdGrr76qKVOmRGy7FCgA\nwLhmmqZKS0u1fPnyYa/4ys/P15IlS2KcbGS9vb3y+XxauHBheJlhGCooKNDhw4dtTDa08Za331NP\nPaVHH300YjcR7jd+ajoAAEN46aWXNHHiRD399NPDrjNjxgylpKTEMNXogsGg+vr6lJSUNGB5UlKS\nOjo6bEo1vPGWV5J27NihtrY2tbS0RHzbFCgAwLjl8/m0fv16tba2jrheXV1dbALBMc6cOaNly5bp\nvffek8fjifj2OYUHABi33n//fV24cEFpaWnyeDzyeDz6+OOPVVFRoVmzZtkdb0SJiYlyu90KBAID\nlgcCASUnJ9uUanjjLa/P59OFCxeUlZUVfm/s379f69at08SJE4edTzdWFCgAwLhVWlqqEydOyO/3\nh3/uvPNOLV++XI2NjXbHG5HH41F2draamprCy0zTVFNTk/Ly8mxMNrTxlregoEAffvih2trawu+N\n++67TyUlJfL7/SNe0TkWnMIDADhad3e3Tp06FT5i0NnZKb/fr4SEBKWlpemOO+4YsL7H41FycrK+\n/vWvh5eVlZUpNTVVq1evjmn20VRUVKi8vFzZ2dnKyclRTU2Nenp6VF5ebne0IY2nvJMmTdI999wz\naNnUqVMjcnsZChQAwNFaWlqUn58vwzBkGIYqKyslXStFQ10+P9SRha6uLrnd7qhntaq4uFjBYFBV\nVVUKBALKzMxUY2Ojpk2bZne0IY23vF90o0ed/isKFADA0ebPn69QKDTm9Ts7Owcta25ujmSkiPJ6\nvfJ6vXbHGLPxlve/iuT7gDlQAAAAFlGgAAAALKJAAQAAWMQcKACAo7S3t0d9jNOnT8dsrEggb/RZ\nzWrc6BdJAQAQCYZhTHe5XB2hUOiWWIzncrksTU63G3mjz+VyfR4Khb5pmuZfR1uXAgUAcAzDMKZL\nSozRcBMlXYnRWJFA3ugLjqU8SRQoAAAAy5hEDgAAYBEFCgAAwCIKFAAAgEUUKAAAAIsoUAAAABZR\noAAAACyiQAEAAFhEgQIAALCIAgUAAGARBQoAAMAiChQAAIBFFCgAAACLKFAAAAAWUaAAAAAsokAB\nAABYRIECAACwiAIFAABgEQUKAADAIgoUAACARRQoAAAAiyhQAAAAFlGgAAAALKJAAQAAWESBAgAA\nsIgCBQAAYBEFCgAAwCIKFAAAgEUUKAAAAIsoUAAAABZRoAAAACyiQAEAAFhEgQIAALCIAgUAAGAR\nBQoAAMAiChQAAIBFFCgAAACLKFAAAAAWUaAAAAAsokABAABYRIECAACwiAIFAABgEQUKAADAIgoU\nAACARRQoAAAAiyhQAAAAFlGgAAAALKJAAQAAWESBAgAAsIgCBQAAYBEFCgAAwCIKFAAAgEX/H/bC\nDFI5JE0LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fd8d67f3860>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from fealpy.functionspace.tools import function_space\n",
    "from fealpy.mesh.TriangleMesh import TriangleMesh\n",
    "\n",
    "#from matplotlib.font_manager import FontProperties\n",
    "#cfont = FontProperties('SimHei')\n",
    "\n",
    "degree = 4\n",
    "\n",
    "point = np.array([\n",
    "    [0,0],\n",
    "    [1,0],\n",
    "    [0,1]], dtype=np.float)\n",
    "cell = np.array([[0, 1, 2]], dtype=np.int)\n",
    "\n",
    "mesh = TriangleMesh(point, cell)\n",
    "V = function_space(mesh, 'Lagrange', degree)\n",
    "ldof = V.number_of_local_dofs()\n",
    "ipoints = V.interpolation_points()\n",
    "cell2dof = V.cell_to_dof()\n",
    "fig, axes = plt.subplots(1, 3)\n",
    "mesh.add_plot(axes[0], cellcolor='w')\n",
    "mesh.find_point(axes[0], showindex=True, fontsize=12, color='g', markersize=25)\n",
    "#axes[0].set_title('local vertices indices', y=1)\n",
    "\n",
    "mesh.add_plot(axes[1], cellcolor='w')\n",
    "mesh.find_point(axes[1], point=ipoints[cell2dof[0]], showindex=True, fontsize=12, markersize=25)\n",
    "#axes[1].set_title(str(degree) + ' dof indices', y=1)\n",
    "\n",
    "axes[2].axis('tight')\n",
    "axes[2].axis('off')\n",
    "rl = [str(i)+\":\" for i in range(ldof)]\n",
    "axes[2].table(cellText=V.cellIdx, rowLabels= rl, colLabels=['m', 'n', 'k'], loc='center')\n",
    "#axes[2].set_title(str(degree) + 'index matrix $I$', y=1)\n",
    "plt.tight_layout(pad=1, w_pad=1, h_pad=1.0)\n",
    "plt.savefig('/home/why/ppt/figures/dof4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 基函数梯度的计算公式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**目标:**\n",
    "\n",
    "给定重心坐标 $(\\lambda_0, \\lambda_1, \\lambda_2)$, 利用**面向数组**的方式计算 $p$ 次元所有基函数在 $(\\lambda_0, \\lambda_1, \\lambda_2)$ 处的梯度值."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**分析:**\n",
    "\n",
    "根据求导法则, 关键是计算出矩阵\n",
    "\n",
    "$$\n",
    "B = \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "1 & 1 & 1 \\\\\n",
    "\\lambda_0 & \\lambda_1 & \\lambda_2\\\\\n",
    "\\prod_{l=0}^{1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\prod_{l=0}^{p-1}(\\lambda_0 - \\frac{l}{p}) & \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p})\n",
    "& \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "中每个元素的梯度. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "构造矩阵:\n",
    "$$\n",
    "F^i = \n",
    "\\begin{pmatrix}\n",
    "1 & \\lambda_i & \\cdots & \\lambda_i \\\\\n",
    "\\lambda_i - \\frac{1}{p} & 1 & \\cdots & \\lambda_i - \\frac{1}{p} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "\\lambda_i - \\frac{p-1}{p} & \\lambda_i - \\frac{p-1}{p} & \\cdots & 1\n",
    "\\end{pmatrix}_{p\\times p}\n",
    ", \\quad i = 0, 1, 2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记\n",
    "$$\n",
    "f_{i,j} = \\sum_{m=0}^j\\prod_{k=0}^j F^i_{k, m},\\quad i = 0, 1, 2, 0 \\leq j \\leq p-1.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\nabla B = &  \n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "\\nabla \\lambda_1 & \\nabla \\lambda_2 & \\nabla \\lambda_3\\\\\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_1 - \\frac{l}{p}) & \n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_2 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{1}(\\lambda_3 - \\frac{l}{p}) \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_1 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_2 - \\frac{l}{p}) &\n",
    "\\nabla \\prod_{l=0}^{p-1}(\\lambda_3 - \\frac{l}{p}) \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0 \\\\\n",
    "\\nabla \\lambda_1 & \\nabla \\lambda_2 & \\nabla \\lambda_3\\\\\n",
    "f_{0,0} \\nabla \\lambda_1 & \n",
    "f_{1,0} \\nabla \\lambda_2 &\n",
    "f_{2,0} \\nabla \\lambda_3 \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} \\nabla \\lambda_1 &\n",
    "f_{1, p-1} \\nabla \\lambda_2 &\n",
    "f_{2, p-1} \\nabla \\lambda_3 \n",
    "\\end{pmatrix}\\\\\n",
    "= & \\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1}  \n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\nabla \\lambda_1 & 0 & 0\\\\\n",
    "0 & \\nabla \\lambda_2 & 0 \\\\\n",
    "0 & 0 & \\nabla \\lambda_3\n",
    "\\end{pmatrix}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "记 \n",
    "$$\n",
    "F =\n",
    "\\mathrm{diag}(P)\n",
    "\\begin{pmatrix}\n",
    "0 & 0 & 0\\\\\n",
    "1 & 1 & 1\\\\\n",
    "f_{0,0} & \n",
    "f_{1,0} &\n",
    "f_{2,0} \\\\\n",
    "\\vdots & \\vdots & \\vdots \\\\\n",
    "f_{0, p-1} &\n",
    "f_{1, p-1} &\n",
    "f_{2, p-1}  \n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用前面的矩阵 $B$, $I$, $Dlambda$, 如下代码即可计算所有单元上的 $p$ 次元基函数的导数:\n",
    "\n",
    "```Python\n",
    "pp = p**p\n",
    "ldof = (p+1)*(p+2)/2\n",
    "Q = B[I,[0,1,2]]\n",
    "M = F[I,[0,1,2]]\n",
    "R = np.zeros((ldof,3), dtype=np.float)\n",
    "R[:,0] = M[:,0]*Q[:,1]*Q[:,2]\n",
    "R[:,1] = Q[:,0]*M[:,1]*Q[:,2]\n",
    "R[:,2] = Q[:,0]*Q[:,1]*M[:,2]\n",
    "gradphi = np.einsum('ij, kj...->ki...', pp*R, Dlambda) # 爱因斯坦求和 \n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<h1><center>谢谢大家！</center></h1>\n",
    "<h1><center>欢迎批评指正！</center></h1>"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
